Interactive comment on “ A numerical scheme to calculate temperature and salinity dependent air-water transfer velocities for any gas ”

i) Intended audience: As pointed out by both referees the scheme put forward does not account for the various effects outside of the standard ’two-phase’ model of gas exchange such as bubble mediated transfer, microlayer effects etc. Necessarily, it also utilizes specific parameterisations for each individual component of the transfer velocity calculation. These are limitations which need to be clearly stated in the abstract and further in the main text of the manuscript, although it is important to consider the uncertainties associated with these paramterisations in the context of the uncertainty in the concentration fields used to drive flux calculations for poorly studied / understood trace gases, which are likely to be large. However, this scheme is not intended to replace either the NOAA COARE algorithm (as both referees clearly recognize) or the detailed study of exchange for a particular gas; rather it is envisaged that the intended users of this scheme will be non-experts in air-sea gas exchange who wish to calculate fluxes of their (probably poorly studied) gas of interest, and where they would currently apply e.g. the Nightingale et al 2000 parameterisation, with a fixed Schmidt number, irrespective of solubilty, temperature etc, and would not even consider the effects of bubbles / chemical enhancement etc. If this scheme is adopted by such studies, then they would at least be employing the standard two-phase model in a consistent, applicable and traceable manner.


Introduction
The rate of exchange of a trace gas between the atmosphere and ocean (or other water surface) is often calculated from observed or inferred concentrations of the gas of interest, from point measurements in individual field studies or averaged or modelled data used in regional or global budgets of trace gas fluxes.When such calculations are undertaken, the "two-phase" model of gas exchange, as applied to the air-sea interface by Liss and Slater (1974) is often used.This model assumes that the transfer of gases is controlled by layers of molecular diffusion on either side of the interface.The flux of a gas across the air-water interface can be expressed as that through either of the two molecular layers (Liss and Slater, 1974) (Eq.1).
where C g and C l are the bulk gas phase and liquid phase concentrations respectively and K H is the dimensionless gasover-liquid form of the Henry's law constant, which is the equilibrium ratio of concentrations in the gas and liquid phase for the gas in question: where C sg and C sl are the inter-facial equilibrium concentrations of the gas in question.K a and K w are the total transfer velocities for the system as expressed from the point of view Published by Copernicus Publications on behalf of the European Geosciences Union.
M. T. Johnson: Air-water gas transfer velocity scheme of the liquid and gas phases respectively.Each are composed of the two single-phase transfer velocities, one representing the rate of water side transfer (k w ) and one representing the air-side transfer (k a ): The transfer velocity in either phase can be represented as an empirical or physically-based function of the physical forcings (typically wind speed, u), the gas-(and medium-) dependent Schmidt number, S c and n, the exponent to which the Schmidt number is raised (representing the degree of renewal of the surface layers due to turbulent mixing), Eq. 5.
It is generally accepted that k w ∝ S −0.5 c and that k a ∝ S −0.67 c (most empirical parameterizations are implicitly "tuned" to these exponents) and these relationships are not investigated here.Various parameterizations of the physical forcing scaling factors for k a and k w are available in the literature (e.g.Liss and Merlivat, 1986;Wanninkhof, 1992;Nightingale et al., 2000).Along with k a and k w parameterizations, three gas-specific terms are required in order to calculate the transfer velocity for a given gas: the Henry's law solubility (which is used in calculating the total transfer velocity from the component velocities), and the air and water Schmidt numbers of the gas.Approaches more detailed than empirical relationships between wind speed and transfer velocity are often used e.g.NOAA COARE (Fairall et al., 1996(Fairall et al., , 2003) ) or MESSY AIR-SEA (Pozzer et al., 2006).However, such generally physically-based schemes require extensive information on the physical properties of the gas and/or more detailed physical forcing data than wind speed alone and are, therefore, more appropriate for well-studied gases such as CO 2 , O 2 or DMS (dimethyl sulfide), where concentration uncertainty is better constrained and physical characteristics of the gases have been determined experimentally, or where a detailed study of physical parameters is commonly undertaken, e.g. in eddy-covariance or other "micro-meteorological" flux studies.Furthermore, such schemes contain tunable parameters which are poorly quantified due to a paucity of data and may vary for different gases due to unquantified or unaccounted for gas-specific effects.While a simpler, broadly empirical scheme cannot hope to address these issues, there is a good case for a generalized scheme which will allow the quick and simple quantification of trace gas fluxes given concentration data, wind speed, temperature and salinity and the minimum physico-chemical information about the gas possible, to fulfil the requirements of a long "tail" of studies covering a wide range of poorly-studied trace gases (from e.g.ammonia to mercury vapour to long-chain hydrocarbons), which are often conducted by biogeochemists or microbiologists rather than experts in quantifying air-water fluxes.
The selection of air-side and water-side transfer velocity parameterizations is the focus of the first section of this paper.Using literature data on solubilities in freshwater and seawater, an empirical relationship between solubility and salinity has been derived (in the absence of a suitable existing relationship) and this is presented in Sect.3. Parameterizations are available for viscosity and diffusivity in air and water and terms from which they are derived, and measured data is available to quantify the Henry's law solubility, and the temperature-dependence thereof, and these are presented in the description of the full numerical scheme in Sect.4, with further details provided in Appendix A. After a description of the software implementation of the scheme (Sect.5), sensitivity analysis and comparison of this scheme with other data/models is presented in Sect.6.The supplementary material contains a full implementation of the scheme, along with input data for 90 trace gases of possible interest to end users, coded in the open-source R software environment, which is freely available for use in other studies.

Selection of transfer velocity parameterizations
Transfer velocity parameterizations are the source of much ongoing investigation and debate in the community and it is not the aim of this work to indicate which one this author believes is "best", or closest to reality (if the two things are different).There is significant uncertainty (commonly quoted as a factor of 2) in the water-side transfer velocity, which is far better studied than k a .It is realistic to assume that the uncertainty in the air-side transfer velocity is, therefore, even greater, with almost complete lack of measurement/validation for trace gases.However, even for dimethylsulfide (DMS), for which there are in the region of 50 000 surface seawater concentration measurements compiled into a recent climatology, uncertainty in the estimated global ocean-atmosphere DMS flux due to the uncertainty in the extrapolated global concentration fields is at least as large as the uncertainty introduced by the difference between using Liss and Merlivat (1986) and Wanninkhof (1992) to estimate the transfer velocity (Lana et al., 2010).Therefore, it can be argued that for pretty much any gas less well-studied than DMS, the selection of transfer velocity parameterizations is probably a small consideration in terms of total uncertainty in understanding fluxes (at least when extrapolating point measurements over large scales of time and/or space).Nonetheless, it is necessary to select a default parameterization for both phases in order to implement the scheme presented here and it is useful to take an informed approach to this selection to aid understanding, if not accuracy.Within the R programme provided in the supplementary information, all of the simple empirical parameterizations mentioned below are implemented along with the physically based model of Woolf (1997); so the "selection" process below is merely a selection of the default parameterization for the scheme (Sect.5 discusses implementation of the scheme and the selection of alternative default sub-models).In the future, it can be expected that improvements to both our physical understanding of gas exchange and our ability to measure trace gas fluxes directly will lead to new and better parameterizations which could and should be implemented in this scheme to supersede those currently present.

Air-side transfer velocity, k a
The air-side transfer velocity is generally expressed as a function of wind speed, u (or specifically, when considering field measurements, u 10 -the wind speed at 10 m above the water surface).A commonly used parameterization in studies of trace gas exchange is that of Duce et al. (1991) [henceforth D91], (Eq.6), which is derived from micrometeorological theory (detailed below) and which depends on the molecular weight of the gas (MW) as well as wind speed.
This is compared in Fig. 1 with the correlations for k a presented in Mackay and Yeun (1983) [MY83], Liss (1973) [L73] and Shahin et al. (2002) [S02] (Eqs.7 to 9) for O 2 and CHI 3 (chosen as examples of relatively small and large molecules, respectively).L73 and MY83 are both derived from wind tunnel studies and S02 from an in situ study using a water-surface-sampling device.Note that there are other wind and field studies of k a resulting in parameterizations which are not considered here, as this analysis is focussed on experimental measurements based specifically on trace gas fluxes rather than water vapour flux, the experimental data which can be considered to be efficiently synthesized by the bulk parameterizations such as D91 and Jeffery et al. (2007Jeffery et al. ( , 2010) ) where S c a is the Schmidt number and D a is the diffusion coefficient of the gas in air (Eqs.23 to 26 in the description of the numerical scheme, Sect.4).The friction velocity, u * is related to the wind speed by the drag coefficient, C D (Eq.10).
MY83 apply the wind speed dependent C D parameterization of Smith (1980) when extrapolating their wind tunnel data to environmental conditions (Eq.11): The wind-tunnel-derived formulation of MY83 overestimates k a relative to D91.As recognised by roughness effects (because of the the non-equilibrium wave field under short fetch) and edge effects (increased turbulence due to wave reflection).Conversely (and in contradiction with most of the literature), S02 suggest that a parameterization based on laboratory measurements should underestimate k a , thus, explaining the fact that their parameterization is approximately 2 times greater than the average of all of the other parameterizations they compare in their paper.This seems an unlikely explanation, and a problem is highlighted when investigating the model they use to validate their data.They find good (1:1) agreement with a resistance model derived from Wesely and Hicks (1977) and Hicks et al. (1987), which are forebears of the parameterization presented by D91 and in most respects very similar.Thus, it is surprising that they should derive a parameterization giving results almost 4 times those of D91 from their observations/model.This seemingly paradoxical result can be explained by the one substantial difference in the aerodynamic resistance term in the S02 model relative to the D91 parameterization: the S02 term is inversely proportional to the variance of the wind direction in radians.The experiments of S02 were conducted on an urban rooftop, where wind direction variability can reasonably be expected to be substantially higher than over a surface with much shorter roughness length and long fetch (such as open water).Typical standard deviation of wind direction in the data presented by S02 was 0.5 radians.A five-fold reduction in this, for instance, brings their parameterization roughly in line with that of D91.The parameterization of Shahin et al. (2002) is, thus, discounted from further consideration as it is probably not applicable to an open water situation.
Further investigation of the D91 parameterization is revealing.The full derivation can be found in their paper, but their parameterization can be summarized in terms of S c a , u * and u 10 (Eqs.12 and 13).
They apply a fixed value drag coefficient, C D (invariant with wind speed) with a value of 1.3x10 −3 , leading to the following relationship between u * and u 10 : Duce et al. (1991) hold C D constant as they consider the wind-speed-dependent component of the drag coefficient to be minor.However, substitution of a wind-dependent term for C D (e.g.Eq. 11) into Eq.12 changes k a substantially at higher wind speeds in the D91 parameterization (Fig. 2).Parameterizations of the drag coefficient are commonly of the form (a + b × u 10 ) 10 −3 Guan and Xie (2003).Smith (1980) is commonly used, but more recent parameter values can be found in the literature, such as that of Yelland and Taylor (1996) (where a = 0.60 and b = 0.07).However, these all tend to be very similar to Smith (1980), causing no significant difference to the calculated transfer velocity.This is not the case with physically derived C D determined from the COARE algorithm, which is discussed below.
The D91 k a term is improved by application of a more realistic term for the drag coefficient, demonstrating a nonlinear response to wind speed, in agreement with the wind tunnel study of MY83.However, it appears to have a much weaker response to changes in the Schmidt number (i.e.difference    2010), demonstrating the sensitivity to the term used for the drag coefficient NOAA COARE drag coefficient data extracted from graph in Fairall et al. (2003).
in diffusivity between different gases) than that suggested by the observations of MY83, who studied a suite of organic molecules covering a wide range of solubilities and Schmidt numbers.
The NOAA COARE algorithm (Fairall et al., 2003) applies a similar resistance model to D91 to calculate k a from wind speed, temperature and atmospheric stability (along with solubility and diffusivity/Schmidt number of the gas in question).The term for k a from the COARE algorithm has been presented specifically for trace gas exchange by Jeffery et al. (2010) [J10] (Eq.14): where κ is the von Karman constant (commonly taken to be 0.4 in seawater).In J10, the conversion from wind speed to friction velocity is via the drag coefficient, but also includes separate terms for gustiness and atmospheric stability, which wind speed / m.s −1 k w_660 / cm.hr −1 q q q q q q q q q q q q q q q qq qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Nightingale et al. ( 2000) Woolf (1997) Wanninkhof (1992)  wind speed / m.s −1 k w_660 / cm.hr −1 q q q q q q q q q q q q q q q q q q q q q q q q q q Fig. 5. Comparison of various observations, parameterizations and tuned model predictions of the increase with wind speed of the liquid phase transfer velocity (k w ) scaled to a Schmidt number of 660.
require various parameter values other than wind speed (inversion height, buoyancy, and β, a tunable parameter) for the calculation.These and other papers presenting modifications of the COARE algorithm specific to gas exchange use the COARE algorithm itself to calculate the drag coefficient and do not present parameterizations of it in their papers (Fairall et al., 2000;Hare et al., 2004;Jeffery et al., 2007Jeffery et al., , 2010)).
Applying one of the conventional parameterizations of the drag coefficient, such as Smith (1980), yields a value of k a with wind speed which is in very good agreement with the D91 equation modified to use the same drag coefficient term (Fig. 3).However, applying the neutral conditions drag coefficient (C D 10N ) presented by Fairall et al. (2003) (taken directly from the COARE algorithm) gives a rather different result (Fig. 4) with substantially lower estimates of k a at low wind speeds (> a factor of 3 difference at 6 m/s).This highlights the potential sensitivity of these k a parameterizations to the drag coefficient and suggests that conventional parameterizations of k a may be substantially overestimating the transfer velocity at low wind speeds, or the COARE algorithm underestimating.However, the conventional terms for the drag coefficient determined from observations e.g.Smith (1980), are correlations between a single parameter (wind speed) and the measured drag coefficient so must implicitly account for other factors, which are represented by discrete, explicit terms in the COARE algorithm.This is clearly a source of uncertainty in the simpler parameterizations but as a simple bulk term it is probably more appropriate to use an empirical parameterization of C D .The J10 parameterization of k a using the Smith (1980) C D term is in good agreement with D91 with the same term for C D applied, but is better at representing the observed spread in k a due to the difference in diffusivity between gases, so it will be adopted as the default parameterization for this scheme with one modification -the addition of the still air diffusive flux observed by Mackay and Yeun (1983) to better represent the zero wind case.The parameterization derived for the numerical scheme presented here is, thus: where u * is calculated from to Eqs. 10 and 11: u * = u 10 6.1 × 10 −4 + 6.3 × 10 −5 u 10 ( 16) The new parameterization presented here is an improvement on previous terms for k a by accounting for the apparent nonlinear relationship with wind speed and the purely diffusive transfer at zero wind and is selected as the default parameterization in the R implementation of the scheme presented in the supplementary material.Note that a relatively recent study has demonstrated that the drag coefficient appears to level out or even start to decrease at hurricane-force wind strengths (greater than ∼40 m s −1 ) (Powell et al., 2003), and this is not accounted for in Eq. ( 16), so further consideration is necessary if calculating fluxes for extreme wind conditions.

Water phase transfer velocity, k W
The water phase transfer velocity (k w ) has been the subject of a much greater study than k a and as such there is an extensive array of parameterizations available, covering a substantial range of possible values and types of relationships. Figure 5 summarises some of the key data and parameterizations of the k w -u 10 relationship as well as two different tunings of the NOAA COARE algorithm.
Parameterizations for k w have been determined empirically in a number of ways, which can be summarised as in situ (deliberate) tracer studies, global (bomb 14 C) tracer studies and direct measurements of trace gas fluxes by micrometeorological (most recently and reliably eddy covariance) methods.All of these are represented in the above diagram as both summary data points and empirical fits.The NOAA COARE algorithm is a fully physically-based model with tunable parameters, the two versions of which shown in Fig. 5 are tuned in different ways.The fit of Hare et al. (2004) is tuned to the Gasex '98 eddy covariance-measured CO 2 fluxes, whereas that of Jeffery et al. (2010) integrates a global ocean modelling approach, bomb 14 C and deliberate tracer data to provide a general tuning, which can be considered the current state-of-the-art for NOAA-COARE.Note the lack of data to validate parameterizations/model tunings at low winds (<2 m/s); and that both of the tuned COARE algorithm lines, along with the cubic fit of McGillis et al. (2001) predict non-zero flux at zero wind speed, whereas all of the other empirical parameterizations (plus the physically based scheme of Woolf, 1997) predict zero transfer velocity at zero wind.There will undoubtedly be a diffusive flux at zero wind through the water surface layer, plus other effects such as rain or convective turbulence may become important; but it is unlikely to be anywhere near that predicted by the fits to higher wind speed data which are not constrained by data at low wind speeds; and any effect is likely to be less significant than in the less viscous air-side surface layer.In the absence of good observational evidence for the value of this diffusive flux, we will assume that a parameterization which passes through the origin is most appropriate.Ultimately, for environmental applications at least, the sensitivity of flux estimates to low wind speed transfer velocities is small because (i) low wind speeds are infrequent, and (ii) fluxes during low wind speeds are small and constitute only a tiny proportion of the totals when integrated over time.
Also included in Fig. 5 is the physically based wind-speeddriven parameterization of k w presented by Woolf (1997) which includes the potential effect of bubbles on the enhancement of transfer, an effect which is likely to be large for insoluble gases and small or non-existent for more soluble gases.Note that the bubble component of this parameterization is the same as that implemented in the NOAA COARE algorithm.There is experimental, and more recently, field evidence for the transfer of more soluble gases being lower at high wind speeds (and greater white-capping), such as for DMS (Huebert et al., 2010) where there is an apparent linear relationship of k w with wind speed (in agreement with a purely turbulent, non-bubble-mediated transfer as predicted by NOAA COARE) and lower predicted k w than Schmidt number scaling of the typical empirical k w -u 10 relationships (which are based on the potentially bubble-enhanced transfer of CO 2 or less soluble gases) would predict.However, this is based on relatively few data points and there are significant uncertainties in bubble-mediated transfer velocity schemes (e.g.Woolf, 1997;Fairall et al., 2003), not least the empirical wind-speed-white-capping relationship, observations of which vary by at least an order of magnitude in the literature (e.g.Moat et al., 2009).Nonetheless, the parameterization of Woolf (1997) is a useful one to constrain the possible effect of bubbles on air-sea gas transfer and is implemented as a non-default option in the R implementation of the scheme.More recent bubble-mediated transfer velocity schemes, which have been validated in the field for noble gases (e.g.Stanley et al., 2009), are not solely dependent on the solubility of the gas in question but also on the partial pressure difference between the atmosphere and underlying waters.The effects of considering disequilibrium between air and water on the significance of bubble effects in trace gas exchange is beyond the scope of the work presented here, but requires further investigation in the future.
The bomb 14 C derived global estimates of k w (e.g.Wanninkhof, 1992;Sweeney et al., 2007;Müller et al., 2008) make the assumption that k w ∝ u 10 2 and scale the k w parameterization to a single point estimated from radiocarbon inventories and globally averaged winds.This analysis is subject to significant uncertainties and as highlighted by Wanninkhof (1992), long-term and short-term averaged winds have different representative mean values leading to different parameter values in the wind speed relationship depending on the application (see Wanninkhof (2009) for the state of the art).In the case of global extrapolations of trace gas fluxes, it may be most appropriate to use a recent bomb 14 C global estimate of k w such as Sweeney et al. (2007) (which applies an improved estimate of the bomb 14 C inventory of the ocean to the Wanninkhof (1992) approach).
The global estimate of k w by Sweeney et al. ( 2007) is in good agreement with the deliberate multiple-tracer study of Nightingale et al. (2000), which is the most comprehensive tracer study published to date.The argument has been made by Nightingale et al. (2000) and others that the dual gaseous tracers used ( 3 He and SF 6 ) may yield a k w that is applicable only to gases of similar solubilities (due to bubble effects) and that such an approach may lead to an overestimation of k w for CO 2 of between 11 and 17% at high wind speeds (assuming the same white-capping-wind-speed relationship employed by the Woolf (1997) parameterization).However, this suggestion is apparently at odds with the good agreement between Nightingale et al. (2000) and Sweeney et al. (2007), highlighting the uncertainty surrounding the magnitude and significance of the bubble effect, at least for CO 2 .For the purposes of this scheme, as the best in situ estimate of k w available in the published literature to date, the parameterization of Nightingale et al. (2000) is selected as the default, with potential 20% uncertainty from the bubble effect being assumed minor relative to the overall uncertainty in calculating the flux of a poorly studied trace gas.Asher (2009) suggests that up to 50% of the variability between different tracer-and wind tunnel-derived k w values can be attributed to experimental uncertainty, further supporting the idea that at the level of our current knowledge of k w the uncertainty introduced by failing to account for the potential effect of bubbles on k w is likely to be relatively small.

Derivation of a salinity-dependent model for Henry's law solubility
In studies of trace gas exchange to date, where measured data on solubility in seawater has not been available for the particular gas of interest, a blanket 20% decrease has commonly been applied to the freshwater solubility to account for the "salting out" effect, following Stumm (1981), irrespective of the properties of the gas in question.The effect of dissolved salts on gas solubility is a subject that remains unresolved, with various competing conceptual models being applied to understand and predict such effects, including ion interaction models (e.g.Pitzer, 1991) and scaled particle theory (e.g.Masterton and Lee, 1970).Whatever model best explains the processes behind the salting out (or sometimes "salting in") of gases, the empirical relationship between solubility and electrolyte concentration can be described according to Setschenow (1889): where ς and ς 0 are the solubilities of the non-electrolyte in the salt solution of interest and in pure water, respectively, C salt is the molar concentration of the salt solution and K s is the empirical Setschenow constant.Here this equation is adapted to represent the salinity-dependent change in solubility as a change in the Henry's law constant: where S is the salinity and K H,0 and K H are the Henry's law constants in pure water and the saline medium, respectively.Note that the negative sign is removed because the gas-overliquid form of the Henry's law constant increases with decreasing solubility.Thus, if K s is known for a particular gas in seawater, it is possible to calculate K H at salinity S, given K H,0 : Extensive experimental data is available for K s for various gases in single electrolyte solutions, but similar data for real or synthetic seawater is considerably more sparse.What data there is can be reproduced reasonably using the numerical scaled-particle theory model of Masterton (1975).However, the use of this model requires gas-specific data on molecular diameter and polarizability, which are not easily available for a wide range of gases and would require additional inputs to the scheme presented here.Therefore, a simpler empirical relationship has been sought.
Rather weak positive relationships between the "LeBas" liquid molar volume of a gas at its boiling point (V b ) (LeBas, 1915) and K s were found in single-electrolyte solutions for a number of gases by Xie et al. (1997) and for a number of gases and liquids by Ni and Yalkowsky (2003).Ni and Yalkowsky (2003) also found a more reasonable linear relationship between K s and the octanol-water partitioning coefficients (K ow ) of the solutes they investigated.
Within scaled particle theory, as applied to the effect of dissolved salts on non-electrolyte solubility, there are two competing processes which can be conceptually related to V b and (K ow ): (i) the salting-out effect of the increasing free energy cost of forming "cavities" in the water matrix to accommodate the non-electrolyte molecules with increasing ionic strength, and (ii) the salting-in effect of the shift in the hydration/dissociation of non-electrolyte molecules in favour of increasing solubility as a result of the increasing polarity of the water-ion matrix with increasing ionic strength (Masterton and Lee, 1970;Masterton, 1975;Zhou and Mopper, 1990;Ni and Yalkowsky, 2003).The former tends to dominate for larger and more hydrophobic molecules and the latter for smaller, more soluble species (Zhou and Mopper, 1990).Molecular size is clearly related to molecular volume, and hydrophobicity/solubility is strongly related to the octanolwater partitioning coefficient.
In the absence of readily available (K ow ) data for all gases, a relationship between K s , Henry's law (i.e.air-water) solubility and molar volume (which is used elsewhere in the transfer velocity scheme) is sought here.A close fit to the K s data for real and synthetic seawater is found by assuming that K s is a function of the natural logarithm of the molar volume (Eq.20) and a solubility-dependent constant of proportionality, θ.
Using laboratory data on K H /K H,0 for seawater/pure water (Table 1) the setschenow constant (K s ) has been calculated according to Eq. ( 18).Then Eq. ( 20) has been solved for θ, using the additive "Schroeder" method of estimating V b from molecular structure (Partington, 1949), which has been shown to have a smaller error than the LeBas method when compared to measured values of V b for a range of compounds (Poling et al., 2001).The Schroeder V b method is used later and is described in the section outlining the complete numerical scheme (see Table 3).We find a strong, nonlinear relationship between θ and ln(K ∅ H ) (Fig. 6), where K ∅ H is the Henry's law constant in pure water at 25 • C. A cubic model is fitted to this with a root mean square error of 8.76×10 −5 (Eq.21).It is logical that K s should be related to the solubility of the gas along with the molar volume: like the octanol-water partitioning coefficient, the Henry's law constant is strongly related to the polarity of the solute and, as observed by Ni and Yalkowsky (2003): "salts produce a continuum that is more polar than pure water ... [and] the more θ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Fig. 6.Regression of the natural logarithm of the Henry's law constant at S = 0 and T = 25 • C with the molar volume scaling factor, θ for the gases listed in Table 1, with the exception of the data of Sing et al. (1999) (marked as solid point on chart) and that of Elliott and Rowland (1993) (see main text for details).
non-polar the solute, the more it is influenced by the polarity of the solvent."21) Table 1 presents predicted K H /K H,0 data for various gases of interest, compared with measured values, including those used to derive the relationship between θ and lnK ∅ H .Note that the deviation from measured values can be attributed both to the error in the derived relationship and also the errors in calculating V b by the Schroeder method, not to mention any possible errors in the experimental measurements.However, agreement is generally good (less than 5% deviation from measured values in all but six cases).This is a satisfactory level of uncertainty relative to previous approaches, particularly considering that the difference between measured values for the same compound from different studies can vary by this amount or more (e.g.CCl 4 and CH 3 Br, Table 1).It is worth noting that the largest deviations are associated with (i) the measurement of NH 3 solubility by Sing et al. (1999), and (ii) the data of Elliott and Rowland (1993).The study by Sing et al. (1999) was conducted at high temperature and extremely high salinity (Note e, Table 1), and is included only to highlight the steep decrease in salting out effect with increasing solubility (in the absence of more appropriate data for soluble gases in addition to the data points for methanal (Zhou and Mopper, 1990) and H 2 O 2 (Bandstra, 2000)).The study by Elliott and Rowland (1993) has been identified by Moore (2000) as underestimating the salting-out effect for CH 3 Cl, and this appears to be supported by the considerably lower value of K H /K H,0 for CH 3 Br compared with that measured by De Bruyn and Saltzman (1997).For the above reasons the data of Sing et al. (1999) and Elliott and Rowland (1993) are not included in the data used to fit the model.It is recommended that where reliable data for the solubility of a gas of interest at the required salinity is available, this data should be considered before applying this predictive equation.

Henry's law solubility
Data on Henry's law coefficients in pure water are readily available for most gases, and where solubility data or transfer velocity calculations are presented in this paper, K H data from the compilation of Henry's law constants presented by Sander (1999) have been used.These data are presented by Sander (1999) as solubilities in pure water at 25 • C, with units of mol L −1 atm −1 .These are converted to dimensionless gas-over-liquid form (Eq. 2), using Eq. ( 22).Interconversions between various other forms of the Henry's law (or Bunsen) coefficients are documented by Sander (1999).
where K H,0 is the dimensionless gas-over-liquid Henry's law constant in freshwater at a given temperature, T (in K); H is the Henry's law constant in mol L −1 atm −1 at 298.15 K, − soln H is the enthalpy change of dissolution, and − soln H R represents the temperature dependence of the solubility (where R is the gas constant).Values of − soln H R are provided by Sander (1999) and other workers.K H,0 can then be scaled for salinity by using Eqs.( 19) to (21).

Air side transfer velocity, k a
The air-side transfer velocity can be calculated using the scheme presented in Eqs. ( 15) and ( 16).Equation ( 15) requires the Schmidt number in air (Eq.23), which is the ratio of the kinematic viscosity of air (υ a ) and the diffusivity of the gas of interest in air (D a ).υ a is in turn the ratio of the dynamic viscosity of air (η a ) and the density of air (ρ a ).
η a and ρ a are calculated according to the scheme of Tsilingiris (2008), applicable to saturated air, which is assumed to be representative of the bottom few mm or less of the atmosphere over a water surface.
where t is the temperature in • C. Values of parameters S V 0 to S V 4 and S D 0 to S D 3 are listed in Table A1.The diffusion coefficients of gases in air (D a ) are calculated according to Fuller et al. (1966): where P is the pressure in atm (assumed to be unity for all calculations presented in this work), V a is the molar volume of air (assumed here to be 20.1 cm 3 mol −1 ) (Tucker and Nelken, 1990).M r is a function of the relative molecular masses of air (M a ), assumed to be 28.97 (Tucker and Nelken, 1990), and of the gas of interest (M b ): Note that this method is defined by Fuller et al. (1966) as being applicable only to pairs of insoluble gases, but it has been successfully applied to the case of insoluble gases in air (Tucker and Nelken, 1990).In a study of binary soluble gas pairs it performed less well than other more complex parameterizations, consistently overestimating diffusion by approximately 23% (Nain and Ferron, 1972).However, as k a is proportional to ( 1 D a ) 2/3 , such error (which is likely to be substantially smaller for a polar gas in mostly non-polar air) has only a small effect on the calculated k a , so it is recommended here that it is applied irrespective of gas solubility.It is worth noting that the sensitivity of the scheme to this term is relatively modest (Table 5).However, the diffusion coefficients predicted using this method for the gases considered here vary over an order of magnitude, e.g., D a (H 2 ) ≈ 0.5 cm 2 s −1 , D a (CHI 3 ) ≈ 0.06 cm 2 s −1 , so it is important to consider this term in the calculation of k a for any particular gas.
Table 4. Input parameter values for test gases used in sensitivity analysis.Note that − H R , the temperature dependence of the solubility, is fixed at an average value of 5000 K −1 1 × 10 3 10 TG5 1 × 10 3 100

Water side transfer velocity, k w
As discussed in Sect.2, the k w parameterization of Nightingale et al. ( 2000) is used by default in the scheme (Eq.28): where S c 600 is a Schmidt number value of 600, which is the typically quoted Schmidt number of CO 2 at 20 • C in freshwater, to which the Schmidt number of the gas of interest under conditions of interest (S c w ) is scaled.Note that Schmidt numbers are very temperature-dependent, so application of this scaling factor is important even for CO 2 .S c w is calculated according to Eq. ( 29): where υ w is the kinematic viscosity of water, D w is the diffusivity of the gas of interest in water, η w is the dynamic viscosity of water and ρ w is the density.The density of saline water, in kg m −3 , is calculated according to Millero and Poisson (1981), according to Eqs. (A1) to (A6) (in Appendix A).

Water viscosity
The dynamic viscosity of water, η w (or η s in the case of water of non-zero salinity) is calculated using the temperatureand salinity-dependent viscosity model/mixing rule scheme of Laliberté (2007a) (Eqs.30 to 32), which requires the compositions of seawater expressed as mass fractions of component solutes, which are derived here from the standard seawater definition of Millero et al. (2008) (Table A2).
The mixing rule follows the form: where η s is the dynamic viscosity of the mixed-solute solution in cP (centipoise; 1 cP = 10 −3 kg m −1 s −1 ), w w is the mass fraction of water in the solution, η w is the dynamic viscosity of pure water in cP.For each solute in the solution, w i is the mass fraction and η i is the dynamic viscosity attributable to the particular solute.Laliberté (2007a) provides www.ocean-sci.net/6/913/2010/Ocean Sci., 6, 913-932, 2010 Table 5. Analysis of the sensitivity of total transfer velocity (K) to changes in the values of key components of the scheme.Note that the estimated uncertainty is intended to represent maximum likely uncertainty in estimation of parameter values from, e.g., model output (for temperature, salinity and wind speed) and from the parameterizations/measurements used for the others.The sensitivity to a change of the same magnitude as the estimated uncertainty is tested.A negative sign denotes a decrease in K in response to an increase in the parameter value.
% Effect on total K of given increase in parameter value Low solubility Medium solubility High solubility Estimated uncertainty a numerically efficient (and for this purpose sufficiently accurate) term for the temperature-dependent dynamic viscosity of pure water (at atmospheric pressure): where t is the temperature in • C. The dynamic viscosity of each component solute, η i is calculated according to Eq. ( 32): where v 1 to v 6 are experimentally-derived empirical constants for each solute.For the common solutes in seawater, values are provided by Laliberté (2007a).Note that Eq. ( 32) is corrected from the original paper according to Laliberté (2007b).Using the ionic composition of standard seawater (Millero et al., 2008) and pairing ions into the common salt constituents found in seawater (Table A2), a temperatureand salinity-dependent viscosity for seawater can then be calculated.This compares well to the fixed-salinity formulation for seawater viscosity of R. C. Hardy at S = 35 (ITTC, 2006) (Table 2).

Diffusion coefficients of gases in water
Diffusivity of a solute in a liquid solvent can be calculated by a number of methods.Using the tabulated data on measured vs. calculated values of diffusivity given by Poling et al. (2001), it can be seen that for the mostly low molecular weight compounds of interest here, and when water is the solvent, the methods of Wilke and Chang (1955) (Eq.33), Hayduk and Minhas (1982) (Eq.34 and 35) and Tyn and Calus (1975) perform best.The method of Tyn and Calus (1975), whilst likely to be the most accurate in the majority of cases (Poling et al., 2001), requires parachor data which is not available for many trace gases so only the other two methods are considered further.The method of Wilke and Chang (1955) [WC55] gives diffusion coefficients in units of cm 2 s −1 : where T is the temperature in Kelvin, M s is the relative molecular mass of the solvent (18.01 in the case of water), η s is the dynamic viscosity of the solvent in cP, is the association factor of the solvent (2.6 for water) and V b is the liquid molar volume of the gas of interest at its normal boiling point (in cm 3 mol −1 ).Hayduk and Minhas (1982) [HM82] present a similar relationship, also in cm 2 s −1 : and The Wilke and Chang (1955) and Hayduk and Minhas (1982) methods are then used to calculate a mean diffusivity for use in Eq. ( 29).

Liquid molar volume at boiling point (V b )
Both methods require the liquid molar volumes (V b ) of the solutes of interest.These can most effectively be calculated from the additive "Schroeder" method (Partington, 1949) (Table 3).This method is generally in good agreement with experimentally-derived values of V b (Poling et al., 2001), but where experimentally-derived values are available for a particular compound it is recommended that these should be used in preference.

Implementation in R
The numerical scheme presented here has been implemented as a program in the open-source R environment (R Development Core Team, 2010), which is freely available on all operating systems/platforms.The program itself ("K calcs Johnson OS.R") is modular, and uses wrapper functions to call each individual parameterization, allowing simple extensibility or inclusion of new parameterizations.For example, the function ka (the air-side transfer velocity) which is called from higher functions or can be called manually (e.g.>ka("NH3",u=6,T=25)) is a wrapper script calling the default selected k a scheme: It can be seen that other existing schemes (which are all included in the code as separate functions which can be called in their own right) can be set as the default by "uncommenting" them and commenting out the currently selected default scheme.

Input files
Two input files are required by K calcs Johnson OS.R: "sw constituent mass fractions.dat"and "compounds.dat".The former reproduces the data presented in Table A2, for calculating the dynamic viscosity of seawater according to Eqs. ( 30) to (32).The latter provides the necessary gasspecific data to drive the scheme.The data required for each gas listed in compounds.dat is detailed in Table 7 6 Sensitivity analysis Table 5 presents the effect of a change in each of the parameters used to calculate the total transfer velocity; the change being of the magnitude of the estimated uncertainty in each parameter.Hypothetical "test gases" are used to cover a range of the two key gas parameters used to drive the scheme: K H and V b .The test gases are defined in Table 4.Note that analysis at different default parameter values (not shown) reveals that the sensitivities change little over wide ranges of temperature (−5 to 35 • C) and salinity (0-35) and only moderately with wind speed, i.e., there are no major nonlinearities in the model that would affect this sensitivity analysis.
The analysis in Table 5 demonstrates that there are different key uncertainties for gases of differing solubilities.For Table 7. Data required for each compound listed in compounds.datincluding an example entry.Note that all except the first three are concerned with calculation of V b .If the V b column has a non-zero value, this value will be used for V b and the structural data will be ignored.soluble gases, key uncertainties arise from (in order of decreasing importance) parameterizations of k a , − H R ,K H and C D .For insoluble gases, key parameters are k w ,D l , η l , V b and ρ l .Assuming some cancelling of errors, it is reasonable to assume that this scheme is accurate to ± 30% for any gas, within the constraints of the simple wind-driven two-phase model and the inherent uncertainty of the parameterizations of the physical forcings of k a and k w (Sect.2).

Application/comparison
The new scheme presented here could be validated against directly measured trace gas fluxes from, e.g., eddy covariance studies, and one such analysis is conducted below.However, looking for direct agreement is to attempt to validate both the scheme and the underlying two-phase model of gas exchange, not to mention the eddy covariance measurements themselves, which is not the aim of this paper.Furthermore, eddy covariance studies have focussed on insoluble trace gases where k w ≈ K w , so it would be validating simply the Nightingale et al. (2000) or other parameterizations available in the scheme.Rather, a demonstrations is presented where the scheme is applied to effectively reproduce the analysis of DMS and CO 2 eddy covariance data and assessment of the bubble effect presented by Huebert et al. (2010) using NOAA COARE.The scheme is also used to recalculate transfer velocities and fluxes from two previous studies which also apply the two-phase model to approach problems of trace gas exchange.wind−speed / m s −1 k w660 or K w(660) / cm hr −1 q q q q q q q q q q q q q q q qq q q qq q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q Heubert et al 2010 eddy covariance DMS derived kw660  Comparison of eddy covariance measurement derived transfer velocities for DMS and CO 2 normalised to a Schmidt number of 660, compared with NOAA COARE and various outputs from the scheme presented here.Note that total transfer velocities (K w ) calculated using this scheme assume a temperature of 15 • C, and do not normalise the air phase Schmidt numbers.q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q −60 −

CO 2 and DMS transfer velocities and the bubble effect
Huebert et al. ( 2010) find a near linear relationship between DMS transfer velocity and wind speed in their eddy covariance measurements of DMS flux.They suggest that this is due to the lack of enhancement of the DMS flux by the bubble effect relative to CO 2 (e.g.Liss and Merlivat, 1986;Woolf, 1997) and demonstrate good agreement between k w inferred from their DMS eddy covariance measurements and the NOAA COARE algorithm (including bubbles) run for DMS (all normalised to a Schmidt number of 660).The scheme presented here can easily calculate Schmidt number normalised water side transfer velocities and these are presented in Fig. 7 for the Nightingale et al. (2000) and Woolf (1997)  The idea of this analysis is not to validate or pass comment on the Huebert et al. (2010) data or findings, simply to demonstrate that reasonably similar results can easily be achieved by application of this scheme rather than NOAA COARE.Comparing the NOAA COARE DMS k w values (solid black line) with those calculated here using the Woolf (1997) k w parameterization (solid light blue) shows reasonable agreement, although with increasing divergence with increasing wind speed.The total transfer velocity calculated using Woolf (1997) for DMS (solid green) is in better agreement with NOAA COARE DMS and is arguably a better fit than the NOAA COARE predictions to the DMS eddy covariance data itself.It also demonstrates the small but important effect of the air-side transfer velocity on reducing the total transfer for DMS.K w calculated in the same way for CO 2 (green dash) agrees well with COARE for CO 2 (black dash), but note that neither of these do particularly well at matching the Hare et al. (2004) eddy covariance measurements at low winds, where Nightingale et al. (2000) seems a better fit.In fact, K w calculated using Nightingale et al. (2000) could be argued to be the best fit through all the eddy covariance data (CO 2 and DMS lumped together).

More soluble (gas-phase controlled) gas exchange
The scheme is used here to calculate Henry's law constants and transfer velocities for ammonia and (along with measured wind speeds and concentrations) to compare the airside controlled ammonia flux with that calculated by Johnson et al. (2008).The results of this analysis are presented in Fig. 8.The main differences between the scheme presented here and that used by Johnson et al. (2008) are that the latter uses (i) the molecular-weight dependent Duce et al. (1991) parameterization of k a , and (ii) a fixed 20% decrease in Henry's law solubility applied to account for the salting out effect.
Note that whilst there is normally good agreement between the methods, there are occasionally up to 25% differences between them (which are associated with high wind speeds and consequent larger fluxes from the nonlinear k a parameterization in the new scheme).There is a systematic bias in the data plotted in Fig. 8 to being above the one-to-one line, which in the negative region of the graph represents greater downward fluxes from the new scheme and in the positive region, smaller upward fluxes.This is due to the salinity correction applied in the new scheme being (at S = 35) approximately an 8% decrease in solubility, compared to the 20% applied by Johnson et al. (2008).

Gases of intermediate solubility
Finally, the relative contribution of the air-side and waterside transfer velocities (k a and k w , respectively) to the total transfer velocity is investigated.The contribution from both sides of the interface to total transfer has been identified as being important for gases of intermediate solubility (Liss and Slater, 1974;Jähne and Haußecker, 1998;Archer et al., 2007).Archer et al. (2007) present data on the percentage reduction in the total transfer velocity (K w ) as a result of including both k a and k w in their calculations of fluxes www.ocean-sci.net/6/913/2010/Ocean Sci., 6, 913-932, 2010 from a marine time-series of iodine-containing halocarbons, compared to the calculation with k w alone.The scheme presented here is used to reproduce these data (Table 6), assuming temperature, wind speed and salinity based on plots of environmental variables in Archer et al. (2007).Agreement is generally very good, which is positive validation for the scheme presented here, particularly considering that Archer et al. (2007) use experimentally determined diffusivities for the gases in question whereas those calculated here are derived from molecular structure.

Conclusions
Presented above is a numerical scheme for calculating the temperature-and salinity-dependent transfer velocity (and component parameters) of any gas across an air-water interface given basic physical chemistry data for the gas (relative molecular mass, K H , − soln H /R), component elements and molecular structure, the ionic composition of the medium and temperature, salinity and wind speed.A fullyimplemented version of the scheme is provided in the Supplementary Material, along with input files for the composition of seawater and a file containing input data for 90 gases of potential interest to workers in the field of trace gas exchange.
The scheme is shown to agree well with previous methods of calculating k a and k w .A novel approach in determining the salting out coefficients which are used to modify the Henry's law coefficients with salinity has improved estimation of this important term significantly over previous approaches.Like most field-based studies of trace gas fluxes calculated from concentration measurements, this scheme does not currently account for processes of chemical enhancement (e.g.Liss and Slater, 1974), or micro-layer effects (e.g.Guitart et al., 2010), which are important in many cases and should be carefully considered by future users of the scheme for their gas/environment of interest.
It is hoped that the scheme presented here and associated software will provide a basis for improved, traceable transfer velocity parameterizations, particularly for less well-studied gases.It is envisaged that users would override default parts of the scheme for their gas of interest where robust measurement data is available to replace parameterized components of the scheme and, furthermore, that key further uncertainties in quantifying gas fluxes would be considered outside of the scheme presented here (e.g., surfactants and other microlayer effects, chemical enhancement, etc.) Table A1.Parameters for calculating the dynamic viscosity (η a ) and density (ρ a ) according to Tsilingiris (2008).These values are for substitution into Eqs.( 24) and ( 25

Fig. 1 .
Fig. 1.Comparison of different parameterizations of k a , the air-side transfer velocity (Eqs.6 to 9) for O 2 and CHI 3 , calculated at 15 • C, where applicable.Note that the parameterization ofLiss (1973) (black line), is neither temperature nor gas specific.For this and subsequent transfer velocity figures (Figs.2 -5) the right-hand panel reproduces the 0-7 m/s wind speed section of the left-hand plot.

Fig. 4 .
Fig. 4. Comparison of the k a scheme presented here with the modified COARE algorithm from Jeffery et al. (2010), demonstrating the sensitivity to the term used for the drag coefficient NOAA COARE drag coefficient data extracted from graph inFairall et al. (2003).

Fig
Fig. 7.Comparison of eddy covariance measurement derived transfer velocities for DMS and CO 2 normalised to a Schmidt number of 660, compared with NOAA COARE and various outputs from the scheme presented here.Note that total transfer velocities (K w ) calculated using this scheme assume a temperature of 15 • C, and do not normalise the air phase Schmidt numbers.

1 Fig. 8 .
Fig. 8.Comparison of the NH 3 ocean-atmosphere fluxes from measured data as calculated by Johnson et al. (2008) and using the scheme presented in this work.
parameterizations along with theHuebert et al. (2010) data and NOAA COARE-calculated DMS and CO 2 bubblemediated transfer velocities, also fromHuebert et al. (2010).Also presented are total transfer velocities (K w ) using S c 660 normalisation for the water side and k a and K H calculated using T = 15 • C. For comparison, theHare et al. (2004) GasEx '98 eddy covariance derived CO 2 transfer velocity values are also plotted.
MY83, laboratory experiments may tend to overestimate transfer velocities relative to environmental conditions, due to fetch and www.ocean-sci.net/6/913/2010/Ocean Sci., 6, 913-932, 2010 Mackay and Yeun (1983)et al. (1991)k a parameterization with that ofJeffery et al. (2010)for CH 3 I and O 2 , both using the Smith (1980) drag coefficient parameterization.The empirical parameterization ofMackay and Yeun (1983)is included for comparison and to demonstrate that, whilst the two physically based parameterizations generally agree very well,Duce et al. (1991)particularly appears to underestimate the degree to which the Schmidt number affects the transfer velocity for different gases.

Table 1 .
K H / K H,0 (salting out factor) for gases measured in natural or synthetic seawater at salinity between 35 and 36 (except where otherwise stated), compared to values predicted by Eqs.(19) to (21).Notes: * -measurement temperature given in parentheses; x -denotes data not used to derive Eq. 21; a -S = 39; b -S = 34, 0.6 • C difference between fresh and salt-water measurements; c -synthetic seawater; d -mean of 0, 15 and 30 • C measurements; e -measurement in S ≈ 400 at 40 • C.

Table 2 .
Comparison of the Laliberté (2007a,b)method for calculating the dynamic viscosity of multiple-solute solutions applied to seawater at S = 35 with the experimentally-derived method of R. C. Hardy(ITTC, 2006)(only valid at S = 35).All viscosities in kg m −1 s −1 .

Table 3 .
Schroeder additive method for calculating V b .For all atoms/structural items a molecule contains, the sum of the increments will give the molar volume.e.g.CH 2 =CH 2 contains 2 carbon atoms, 4 hydrogen atoms and 1 double bond so the Schroeder V b is 2×7+4×7+7 = 49cm 3 mol −1 * applies to all kinds of cyclic features and is applied only once to ring-containing compounds irrespective of the number of rings present.

Table 6 .
Archer et al. (2007)duction in calculated total transfer velocity (K w ) for a series of iodocarbons from the dataset ofArcher et al. (2007)when considering both k w and k a , relative to the transfer velocity calculated using k w alone.

Table A2 .
). Parameters for calculating solute viscosity, and mass fraction of each solute (as a proportion of total salinity) for Eqs.(30) and (32).