Numerical tools to estimate the flux of a gas across the air-water interface and assess the heterogeneity of its forcing functions

abstr A numerical tool was developed for the estimation of gas fluxes across the air water interface. The primary objective is to use it to estimate CO2 fluxes. Nevertheless application to other gases is easily accomplished by changing the values of the parameters related to the physical properties of the gases. A user friendly software was developed allowing to build upon a standard kernel a custom made gas flux model with the preferred parametrizations. These include single or double layer models; several numerical schemes for the effects of wind in the air-side and water-side transfer velocities; the effect of turbulence from current drag with the bottom; and the effects on solubility of water temperature, salinity, air temperature and pressure. It was also developed an analysis which decomposes the difference between the fluxes in a reference situation and in alternative situations into its several forcing functions. This analysis relies on the Taylor expansion of the gas flux model, requiring the numerical estimation of partial derivatives by a multivariate version of the collocation polynomial. Both the flux model and the difference decomposition analysis were tested with data taken from surveys done in the lagoonary system of Ria Formosa, south Portugal, in which the CO2 fluxes were estimated using the IRGA and floating chamber method whereas the CO2 concentrations were estimated using the IRGA and degasification chamber. Observations and estimations show a remarkable


Introduction
The appropriate algorithms for the estimation of gas fluxes across the air-water interface have been the subject of great concern by the scientific community.One of its most notorious applications is in studies about the CO 2 exchange between the atmosphere and the global oceans (Takahashi et al., 2002(Takahashi et al., , 2009) ) coastal oceans (Frankignoulle, 1988;Frankignoulle and Borges, 2001;Sweeney, 2003;Vandemark et al., 2011), estuaries (Carini et al., 1996;Raymond et al., 2000;Hunt et al., 2011;Oliveira, 2011;Oliveira et al., 2012), rivers (Cole and Caraco, 2001), lagoons (Thomaz et al., 2001) and lakes (Cole and Caraco, 1998;Koné et al., 2009;Sahlée et al., 2011).The marine and aquatic environments may work as either net sinks or net sources of CO 2 for the atmosphere.Nevertheless, this shows a great spatial and temporal variability (Smith and Hollibaugh, 1993;Duarte and Prairie, 2005;Borges, 2005;Borges et al., 2005).The flux of CO 2 across the air-water interface is fundamental to estimate the carbon budget of marine and aquatic ecosystems and classify them as either autotrophic, upon net CO 2 consumption by primary producers, or heterotrophic, upon net CO 2 production by bacterial degradation of organic carbon.Coastal oceans and riverine systems are believed to be globally heterotrophic, remineralizing organic carbon imported from terrestrial ecosystems (Smith and Hollibaugh, 1993;Cole and Caraco, 2001;Duarte and Prairie, 2005;Borges et al., 2005).Although occupying a small fraction of the global ocean, the coastal oceans are major sources of CO 2 to the Published by Copernicus Publications on behalf of the European Geosciences Union.
V. M. N. C. S. Vieira et al.: Air-water interface gas flux atmosphere, presenting an average CO 2 flux per unit area about 5 times higher than the open ocean (Smith and Hollibaugh, 1993).The flux of a gas across the air-water interface has also been studied for the cases of volatile pollutants, such as organochlorine pesticides, hydrocarbons and heavy metals.These are often imported from industrial and agricultural catchment areas through river basins to the coastal waters of highly populated coastal areas where they may be released to the atmosphere.
There are many physical, chemical and even biological aspects mediating the fluxes of gases across the air-water interface.There is also extensive literature covering the majority of these aspects.However, very few attempts have been made to try and integrate several of these factors, particularly when it involves combining distinct fields of knowledge such as chemistry, physical oceanography, meteorology and numerical modelling.Therefore, the first objective of the current work was to develop a numerical tool that provides an accurate estimate of the flux of a gas across the air-water interface.Focus was kept on CO 2 .For studies of other gases, substituting the adequate parameters in the model is required.This numerical tool was based on that of Johnson (2010) but underwent several upgrades: (i) it is possible to choose between single or double layer models; (ii) new numerical schemes for the effect of wind in the water-phase transfer velocity by Mackay and Yeun (1983), Carini et al. (1996), Raymond and Cole (2001), Zhao et al. (2003) and Borges et al. (2004b) were introduced; (iii) the effect of sea surface agitation in the water-phase transfer velocity was added; (v) the effect of atmospheric stability in both the air-side and water-side transfer velocities was added; and (iv) the effect in the water-phase transfer velocity of turbulence due to current drag with the bottom following O' Connor and Dobbins (1958) was added.The latter may play a fundamental role in regulating the gas transfer velocity in macro and mesotidal estuarine and lagoonary systems.The second objective of this work was to develop a numerical method that allows decomposing a difference in the gas fluxes between two distinct situations into the effects of their differences in the environmental variables.This enables the identification of the variables responsible by differences in fluxes between two situations.The current work is intended to set the grounds for further research.This shall consist of including more environmental processes, improving the algorithms of the currently included ones, submitting the tools to a wide range of environmental conditions and conjugating them with numerical modelling labs such as MOHID (www.MOHID.com),ECO lab (www.dhisoftware.com)and WASP (www.epa.gov/athens/research/wasp.html).

State of the art
The flux (mol m −2 s −1 ) of a gas across the air-water interface is usually estimated as F = k(C a /k H − C w ), where k (m s −1 ) is the transfer velocity which often has incorporated the chemical enhancement factor α (scalar), C a and C w (mol m −3 ) the CO 2 concentrations in the air and water, respectively, and k H (scalar) Henry's constant in its C a /C w form.Here, a positive F represents a flux from the air to the water.The gas flux is frequently estimated by the alternative formulation F = kα pCO 2 , where the CO 2 concentration in the air is given in its partial pressure and the CO 2 concentration in the water is given in its expected air partial pressure would it be at equilibrium with the water and α is Bunsen's gas solubility coefficient, equivalent to Henry's constant (k H ) in its C w /P a form.Sander (1999) and Johnson (2010) proposed algorithms to estimate Henry's constant and convert it into its several forms.In order to estimate the effects of water temperature, salinity, air temperature and pressure on solubility/volatility, these formulations consider the physical and molecular properties of the air, gas, water and its solutes.
The k term represents the transfer velocity (also known as piston velocity) of the gas molecules across the air-water interface.In still air and still water conditions, this movement of molecules across the thin layer is due to diffusive transport and thus constrained by the environmental variables that regulate diffusivity.However, when at least one of the phases is not still, turbulence at the interface becomes the main factor regulating the gas transport.The simpler models for the estimation of the transfer velocity consider a single thin layer (Carini et al., 1996;Raymond and Cole, 2001;Borges et al., 2004b;Zappa et al., 2007) across which the transfer velocity equals the water-phase transfer velocity (k = k w ).Full explanation of all the algorithms for the k w estimates would be too extensive and beyond the scope of the current work.In this work, focus is kept in the fundamental physical aspects and the methods to simulate them.A provisional turbulencedriven water transfer velocity (k # w ) is usually estimated as a function of the wind speed (u 10 ) at 10 m height or, alternatively, of the air-side friction velocity (u * ) at the air-water interface (Mackay and Yeun, 1983;Zhao et al., 2003).Most often, these are first to second degree polynomials.A constant with the value of 10 −3 is sometimes added to the k # w representing the transfer velocity in still conditions, i.e. the transfer velocity due to diffusivity when wind speed is zero.There are more physical phenomena that affect the water-side transfer velocity and for which there have been proposed algorithms to simulate them.Such are the cases of the formation of bubbles with high wind speeds and breaking waves (Memery and Merlivat, 1985;Woolf, 1997Woolf, , 2005;;Zhao et al., 2003;Duan and Martin, 2007), wave field (Taylor and Yelland, 2001;Oost et al., 2002;Fairall et al., 2003;Hwang, 2005;Zhao and Xie, 2010), rain (Ho et al., 2004;Zappa et al., 2009;Turk et al., 2010), surfactants (Frew et al., 2004) and the variability of the wind velocity over longer time intervals (Wanninkhof, 1992).The parameterization by Fairall et al. (2000) attempts to congregate the fundamental environmental factors over the open ocean.
The provisional water transfer velocity (k wind w ) is estimated for fresh water at 20 • C and rectified to the final water transfer velocity (k w ) at actual temperature and salinity multiplying it by the chemical enhancement factor (α).This factor is usually taken as (Sc w /600) −0.5 : where Sc w is the Schmidt number of water estimated for the actual temperature and salinity, with 600 usually accepted as the Schmidt number for fresh water at 20 • C and distinct exponents have been proposed, particularly when related to sea surface agitation or the presence of surfactants.The Schmidt number at actual water temperature and salinity may be given by algorithms of a statistic nature (Carini et al., 1996;Raymond and Cole, 2001;Borges et al., 2004b).These are polynomials that best fitted observations.Alternatively, Johnson (2010) proposed a mechanistic numerical scheme that accounts for the effects of temperature and salinity considering several the physical properties of pure water, its solutes, and the diffusing gas.In such a case the mass diffusivity in the water may be estimated by the algorithms proposed by Hayduk and Laudie (1974), Hayduk and Minhas (1982) and Wilke and Chang (1955).Borges et al. (2004b) proposed adding to the wind-driven turbulence the turbulence due to the water current and its drag with the bottom (k current w ) as this may be an important source of turbulence in coastal waters.Its algorithm is given by O' Connor and Dobbins (1958).Woolf (2005) further proposed splitting the k wind w term into a term for sea surface agitation plus a term for whitecap (i.e.bubble formation from breaking waves).
Equation ( 1) is one of the most used formulations for the water-side transfer velocity.It was the adopted in this work and thus was presented with detail.There are nevertheless other two widely used formulations.The bulk model was implemented in the COARE algorithm (Fairall et al., 1996;Grachev and Fairall, 1997;Fairall et al., 2003) to estimate the fluxes of heat, humidity and gases across the air-water interface, forced by wind, atmospheric stability and seasurface agitation, and associated with the eddy-covariance field methodology.Surface renewal theory and micro-scale wave breaking congregate a vast body of literature, developed by B. Jähne, E. J. Bock, and associates at the University of Heidelberg and C. J. Zappa, N. M. Frew, W. R. McGillis and associates at the Woods Hole Oceanographic Institution, devoted to the estimation of the transfer velocities of gases, heat and humidity sustained on a common numerical scheme.The work by Frew et al. (2004), relying on such a scheme, bonds the effects of the main related environmental factors.
A slightly more complex model, the thin film model (Liss and Slater, 1974;Johnson, 2010), also called the tworesistance model (Mackay and Yeun, 1983), considers along the air-water interface both the water-phase and the air-phase thin layers.The final transfer velocity is the weighted har-monic mean of the air-side and water-side transfer velocities, k a and k w respectively.Depending on whether the flux is being estimated from the air-side or the water-side point of view, the transfer velocity scheme weights the opposite phase transfer velocity by Henry's constant.In the formulation above the flux (F ) is estimated from the water point of view and thus the transfer velocity (k) is estimated as in Eq. ( 2).
To compute the flux from the air point of view , the transfer velocity must also be estimated from the air point of view in Eq. ( 3).Despite the different transfer velocities, the fluxes yielded by both methods are equal.
In this thin film model the water-phase transfer velocity (k w ) is estimated, likewise the transfer velocity in the single thin layer model in Eq. ( 1), whereas the air-phase transfer velocity (k a ) needs a different formulation.The k a is mainly driven by the wind velocity.Therefore, Duce et al. (1991), Liss (1973) and Shahin et al. (2002) estimate k a directly from u 10 , whereas Mackay and Yeun (1983), Zhao et al. (2003) and Johnson (2010) estimate it from the friction velocity (u * ).
The simplest way to get to u * from u 10 is through the drag coefficient: CD = (u * /u 10N ) 2 , where N stands for neutral atmospheric stability conditions.The simplest formulation is by Duce et al. (1991) proposing a fixed value drag coefficient, which has been proved to be unrealistic.A variable drag coefficient dependent on u 10 was estimated from field surveys (Smith, 1980), wind tunnel experiments (Mackay and Yeun, 1983) and deep water wind seas (Taylor and Yelland, 2001).Sethuraman and Raynor (1975) proposed drag coefficients dependent on the surface roughness and estimated by the Reynolds number, or dependent on the atmospheric stability and estimated by the Richardson number.Air temperature and pressure may also affect the air transfer velocity.Therefore, Mackay and Yeun (1983), Shahin et al. (2002) and Johnson (2010) propose air transfer velocity equations that include temperature-and/or pressure-dependent terms of the air diffusivity (D a ) and/or the Schmidt number of air (Sc a ).

The gas flux model
The current work provided a numerical scheme for the estimation of the flux of a gas through the air-water interface, a Matlab ® -based free open source software package to implement it and a tutorial for the software (available in the Supplement).Model implementation followed the section above.

V. M. N. C. S. Vieira et al.: Air-water interface gas flux
It could be either as single layer or double layer (thin film), and the transfer velocity estimates relied on Eq. (1).A thorough explanation on the available options is presented in the software tutorial.Several of the k w and k a algorithms relied on the friction velocity, which could be estimated from u 10 using the CD.The most comprehensive model implementation included the effects of sea-surface roughness and atmospheric stability on the turbulence-driven transfer velocities.Surface roughness is dependent on the wave field, and therefore on the wind intensity and on the distance it has been acting upon the water surface (i.e. the fetch) generating a shear stress.The formulation proposed followed the same rationale as the AERMOD, developed by the Environmental Protection Agency (EPA).The basic principle was to adapt the log-wind profile equation solving for friction velocity as a function of wind speed and roughness length.Then, apply this estimate to the available friction velocity-based formulations of air-side and water-side transfer velocities, k a and k w , respectively.Still, atmospheric stability may also play an important role in the relation of wind speed with friction velocity.Thus, a more accurate formulation is the log-linear wind profile in Eq. ( 4), named so because it incorporates a logarithmic term for roughness length and a linear term for atmosphere stability.
Here, u z (m s −1 ) is the wind velocity at height z (m), u s (m s −1 ) the collinear component of the water current velocity at the sea surface, k the von Kármán constant (usually 0.4) and z 0 (m) the roughness length.To avoid confusion it must be noted that z is height in meteorology, whereas depth in oceanography and that von Kármán constant (k) should not be confounded with the transfer velocity (k).The linear term is the atmospheric stability function ψ u , also called ψ m as it is the vertical transfer of momentum being addressed.There are also ψ h and ψ q for the vertical transfers of heat and humidity.It is the integrated non-dimensional gradient given by for atmospheric stable conditions, whereas for unstable conditions are proposed several different algorithms, all much more complex (Businger et al., 1971;Dyer, 1974;Stull, 1988;Lee , 1997).The β is a constant usually between 4.5 and 5, and L (m) is the Monin-Obukhov length given by where u * is the friction velocity (m s −1 ), ρ the air density (g m −3 ), the potential temperature of air (K), c p the specific heat of air (J g −1 K −1 ), H the vertical heat flux (J m −2 s −1 ) assumed positive upwards and g the gravitational acceleration (m s −2 ).Given the absence of data about the vertical heat flux, L was estimated from the bulk Richardson number following two different formulations by Stull (1988) and Lee (1997).This required data at two different heights.It was arbitrarily chosen that at z = 10 air temperature equalled T a , whereas at z = 0 air temperature equalled the temperature of the sea surface.However, the "cool skin" and "warm layer" formulations were not implemented.Air pressure was always P , and relative humidity was 0.7 at z = 10 and 1 at z = 0, following Lange et al. (2004).Potential and virtual temperatures were estimated according to Stull (1988).The log-linear wind profile in Eq. ( 4) was solved for the friction velocity: It was also implemented an alternative formulation by Sethuraman and Raynor (1975) estimating the drag coefficient as a linear function of the difference between air and water temperatures.However, it was of limited use as (i) it is only valid for −3 • C ≤ T ≤ 3.7 • C and (ii) it cannot be conjugated with roughness length formulations.
The field estimates of roughness length were done according to Taylor and Yelland (2001) formulation, z 0 /H s = A(H s / L p ) B , where H s (m) is the significant wave height, L p (m) the wave length of waves at the peak wave spectrum, a scaling constant presently introduced, A = 1200 and B = 4.5.This parameterization predicts the drag coefficient (and thus also the friction velocity and roughness length) increases with increasing fetch and wind duration.Other parameterizations by Donelan (1982Donelan ( , 1990)), Smith et al. (1992), Oost et al. (2002) and Fairall et al. (2003) estimate the wave age based on peak wave speed and friction velocity.These were not tested as their requirement for a friction velocity input would return a circular function.

Field estimates and unit conversions
The wave field data were collected by Instituto Hidrográfico's buoy located 6.1 km off shore from Ria Formosa and over 93 m depth.Boat trips were performed inside Ria Formosa and at the nearby coastal ocean to collect the remaining data.An Excel worksheet with the data is provided together with the software in the Supplement.This worksheet is proposed as the protocol for the required data.
The gas concentrations are commonly estimated from the field in either mol m −3 or ppm units.In the current work were used data with the infrared gas analyzer (IRGA) and floating chamber sampling procedure, yielding the gas concentrations in ppm.The software accepted gas concentrations in either form and converted these into mol m −3 to estimate the fluxes.There were two distinct types of conversions: (i) the [gas] in the air converted between ppm and mol m −3 using the ideal gas law, and (ii) the [gas] in the water converted between mol m −3 and its equivalent air ppm at equilibrium, using Henry's constants.The details on these conversions are provided in Supplement A, together with the protocol for the estimation of the flux from the floating chamber data.Preliminary tests with the model yielded a flux even when the CO 2 concentrations (both given in ppm) in the water and in the air were in equilibrium.It enlightened the need for careful, accurate conversion between the distinct forms of Henry's constants.The k Hpc is Henry's constant for water at 25 • C and 0 ppt salinity given in its P a /C w form.It has a value of 29.4118.The k H is Henry's constant for a given temperature and salinity in its C a /C w form.Johnson (2010) presents an algorithm to estimate k H from k Hpc .This algorithm is represented in the first line of the braced expression in Eq. ( 8).The f (T ) and f (S) represent the functions that resolve for the given temperature and salinity respectively, T K,w the water temperature in kelvin and α H a constant with the value of 12.2.This constant is given by Sander (1999) in an algorithm to estimate k Hcp from k H (in the second line of the braced expression in equation 8) where T K,a is air temperature in kelvin.The k Hcp is needed to convert the equilibrium CO 2 concentration in the water from ppm to mol m −3 at the given environmental conditions (Eq.A5).However, it is fundamental that the k Hcp estimation for those environmental conditions follows the same algorithm previously used for the k H estimation for the same environmental conditions in Eq. ( 8).Furthermore, it is also essential to note that the temperature in Sander (1999) expression is relative to air.This is not explicit in the original article, and one may easily be misled assuming it is water temperature, because this is the main control of solubility.However, its effect was already accounted for in the k H estimation from k Hpc .This is demonstrated by developing the flux equation to C a /(C w • k H ) = 1.If both CO 2 concentrations are given in ppm and their conversions are introduced into this equation, knowing that P (atm) = 101 325.01 Pa, in Eq. ( 9) is obtained; but only if the temperature in Sander (1999) expression is air temperature.Otherwise, the equation only applies when air and water temperatures are equal.Equation (9) was also used to accurately determine α H as 12.1866.

Decomposition of the difference in the gas fluxes (DDF)
For some studies it may be useful to compare a particular case of a gas flux with that of a reference situation, identifying and ranking the causes for the difference.The environmental conditions of the reference situation were recorded in a column vector x a , and its CO 2 flux was estimated by the numerical model above as f a .The environmental conditions of the particular case were recorded in a column vector x b , and its CO 2 flux was estimated by the numerical model as f b .The difference between the environmental conditions of the particular case and of the reference situation ( x) was given in the column vector h in Eq. ( 10).It was chosen to separate atmospheric stability from the remaining effects of air and water temperatures (the software allows doing so).Therefore, the column vectors were arranged as x 1 = C air , x 2 = T air , x 3 = P , x 4 = u 10 , x 5 = z 0 , x 6 = ψ u , x 7 = C w , x 8 = T w , x 9 = S, x 10 = w and x 11 = z w .It is important to note subscript a presently stands for the reference situation and not for air.
The difference in the CO 2 flux was given by f b − f a .It was decomposed into its multiple parcels, each attributable to the difference in a particular environmental variable or interactions between variables.This decomposition was possible developing the Taylor expansion of the gas flux model: There were two sources of error here.One was the difference between observed (f obs ) and estimated (f est ) fluxes.This was only addressed by the gas flux numerical model and not by the DDF.The other was the difference between the estimated left-hand side and the estimated right-hand side of Eq. 11.
This was the remainder of the Taylor expansion tending to zero as tends to ∞.The integer stated the highest order terms used, usually high enough for the remainder to be close to zero.However, as there were many independent variables, the number of higher order terms was too big and its estimation computationally too heavy.Therefore, the software enabled automatically adjusting this decomposition for a specified number of independent variables, each with its own iorder terms in Eq. ( 12).Each term of the Taylor expansion was located in a specified entry of a data array (named Taylor Array) with i dimensions in Eq. ( 13).In this case it was a hypervolume with 11 dimensions.The coordinate of each term in each dimension was given by the respective rank of its partial derivative.This procedure enabled a variable-wise sorting out of insignificant terms, optimizing computational effort.The multivariate form of the Taylor expansion has each term preceded by a coefficient given by the multinomial in Eq. ( 14).However, the numerator in Eq. ( 14) cancels out with the denominator from the middle quotient in Eq. ( 13), thus simplifying the calculus of Eq. ( 15).Subtracting f a was done setting the first entry in Taylor Array to zero.

+1
n 11 =1 Taylor Array n 1 ,n 2 ,...,n 11 (15) The partial derivatives were estimated numerically at point k located within the interval between x a and x b .While the detailed explanation on the procedure is available in Supplement B, here only a brief overview is presented.The gas flux function was approximated by a collocation polynomial in its turn estimated by a multivariate adaptation of Newton's finite difference formula.The collocation polynomial was partially derived to each of the dimensions.The output was a numerical estimate of the partial derivatives of the collocation polynomial that fitted with accuracy the partial derivatives of the gas flux function for any particular point in the hypervolume of independent variables.
Ideally, the whole gas flux difference was partitioned between the independent variables and not between combinations of these variables.To achieve this, each multivariate term of the Taylor expansion was itself evenly partitioned among the independent variables contributing to it.The remainder was estimated by subtracting the sum of the estimated terms from the actual gas flux difference given by f b −f a .It allowed tracking the accuracy of the results, which was one of the criteria used for model optimization.The other was the computational time required to perform the calculus.The model optimization was tested for each dimension at a time and included three features: 3.In the process of numerically estimating derivatives, it is crucial the size of the steps taken forward or backward (the δ i ) in Newton's finite difference formula for the collocation polynomial.If these are too large or too small, with increasing order of the terms, the δ i raised to higher powers leads towards infinity or infinitesimal, which turns the error unbearable.A simple, direct answer to this problem was choosing the δ i to always be in the vicinity of 1.However, for some variables, their increase in steps of size 1 would get them out of bounds, that is, far out of the interval given by x i,a and x i,b .Thus, it was also necessary to play with the units upon which the steps were taken so that they would be within bounds but still represented by numbers with one digit: (a) gas concentrations could be converted from mol m −3 into mmol m −3 ; (b) air pressure from atm into kilopascal (KPa); (c) wind speed from m s −1 into Km h −1 ; (d) roughness length from m into dm, cm, mm or 10 −1 mm; (e) ψ u , a scalar, into •10, •100 or •1000 units, (f) current speed from m s −1 into dm s −1 , cm s −1 , m min −1 , hm h −1 or km h −1 ; and (g) depth from m into dm, dam or hm.
This analysis presented a bias when the [gas] was supplied in units of ppm to a model that works with units of mass volume −1 , and there was a temperature and/or pressure difference between reference and alternative sites.To clearly illustrate this issue, consider a reference and alternative sites that were equal in every variable except air pressure.In this case the reference and alternative sites have equal [gas] when expressed in units of ppm but different [gas] when expressed in units of mass volume −1 , simply because equal amounts of gaseous mass occupy different volumes when subject to different pressures.The bias was not being considered the effect on the gas flux of this [gas] difference induced by the air pressure.Therefore, there was a part of the flux that was missing.Therefore, the numerical estimates of the partial derivatives had to be rectified: when the [gas] was given in ppm, it was not automatically converted to mol m −3 .First, the steps further were taken in Newton's finite difference formula with the [gas] still in ppm units as these were equally well suited for that purpose.Only after each step was taken, the respective ppm was converted to the mol m −3 that was fed to the flux model.This procedure enabled accounting for the effects of air temperature and pressure variations on the conversion of the gas concentrations.

Air-side transfer velocity
Wind was the most influential factor affecting the air-side transfer velocity.Several algorithms simulating this relation are presented in Fig. 1.All the equations about the wind effect including a term for the drag coefficient (Johnson, 2010;Mackay and Yeun, 1983) were very coherent among each other.As expected, the Duce et al. (1991) constant drag coefficient underestimated the air transfer velocity at high wind speeds relative to the drag coefficient parameterizations by Smith (1980) and Mackay and Yeun (1983).Furthermore, this parameterization passed through the origin, meaning no CO 2 flux at still air.Other formulations presented the same problem, as was the case of the COARE formulation by Jeffrey et al. ( 2010).In the COARE algorithm this was solved with the addition of a gustiness term (Grachev and Fairall, 1997;Fairall et al., 2003).Presently, this was solved with the addition of a constant (10 −3 ) following Mackay and Yeun (1983) and Johnson (2010).Only the algorithm by Shahin et al. (2002) simulated perceptible effects of air temperature and pressure on k a (not shown), but these were relatively meaningless.Some formulations use the friction velocity instead of u 10 .This subject is explored in the "friction velocity" section.

Water-side transfer velocity
The water diffusivity equations yielded approximate results with water temperature changing from 0 • C to 30 • C (Fig. 2a).Thus, choosing different diffusivity equations had little effect on either the Schmidt number of water (Sc w ) or the chemical enhancement factor (α) when estimated according to Johnson (2010).Other α algorithms by Borges et al. (2004b), Carini et al. (1996) and Raymond and Cole (2001) also yielded approximate results with changing temperature (Fig. 3a).The estimates of the effect of salinity in both the water diffusivity (Fig. 2b) and α (Fig. 3b) were very subtle.The Carini et al. (1996) and Raymond and Cole (2001) algorithms do not account for salinity in the α estimate.When comparing the several available algorithms for the relation of wind speed with k wind w , two groups were set aside (Fig. 4).The first group had the algorithms developed for open ocean estimates and/or strong winds.Their relations were exponential.The formulation by Mackay and Yeun (1983) was estimated in a wind tunnel with wind speeds between 5 and 22 m s −1 and extrapolated for environmental conditions using the wind-dependent drag coefficient scheme by Smith (1980).The second group had the algorithms developed from river and estuarine surveys in low wind regimes.The Carini et al. (1996) and Borges et al. (2004b) functions were linear.The Raymond and Cole (2001) function is an exponential function estimated exclusively from wind speeds below 8 m s −1 .Its extrapolation to high winds was a wild guess yielding the fastest transfer velocities.Some formulations use the friction velocity instead of u 10 .This subject is explored in the "friction velocity" section.Only one algorithm, by O' Connor and Dobbins (1958), was used to estimate the effect of water current and depth on the water-side transfer velocity (k current w ).The transfer velocity increased nonlinearly with increasing water current and decreasing depth (Fig. 5).Its magnitude was similar to the magnitude of the water transfer velocity imposed by low to moderate winds."H&L74": Hayduk and Laudie (1974); "H&M82": Hayduk and Minhas (1982); "W&C55": Wilke and Chang (1955).

Friction velocity
Some formulations for both the air-side and water-side transfer velocities rely on u * rather than on u 10 , thus allowing to account for the effects of roughness length and atmospheric stability.In Fig. 6 T a was simulated changing from 0 • C to 40 • C while T w was fixed at 15 • C. It simulated atmospherically unstable conditions when T a < T w and atmospherically stable conditions otherwise.The bulk Richardson number was negative in the former case and positive in the latter.The Stull (1988) and Lee (1997) formulations largely mismatched.However, both predict (i) values within the expected bounds, (ii) friction velocity increasing with conditions changing from stable to unstable (increased momentum transfer across the atmosphere boundary layer) and (iii) friction velocity increasing with roughness length (higher roughness lengths created more wind drag).The estimates with the Sethuraman and Raynor (1975) formulation were Fig. 3. Effects of T w (a) and S (b) on the chemical enhancement factor (α) in fresh w Borges et al (2004); 'Cea96': Carini et al (1996); 'R&C01' : Raymond and Cole 2001.J α estimate with water diffusivity by: 'H&L74': Hayduk and Laudie (1974); 'H&M82' Minhas (1982); 'W&C55': Wilke and Chang (1955).Carini et al. (1996); "R&C01": Raymond and Cole (2001), Johnson (2010) α estimate with water diffusivity by: "H&L74": Hayduk and Laudie (1974); "H&M82": Hayduk and Minhas (1982); "W&C55": Wilke and Chang (1955).
only plotted for T within its specific bounds.In this case it corresponded to 12 • C ≤ T a ≤ 18 • C.After wind, roughness length and atmospheric stability were the next most influential factors in both the air-side and water-side transfer velocities (Fig. 7).The schemes by Mackay and Yeun (1983) were fit to wind tunnel data and, thus, in the absence of long fetches (and therefore of rough surfaces) and under neutral atmospheric conditions.However, when these effects were added, the transfer velocity estimates increased significantly and could even surpass the highest estimates by formulations based on oceanic data.

The gas
The overall influence of water temperature and salinity on the CO 2 flux was estimated with the temperature set to 17.  (1985).Where applicable the u * was estimated from u 10 using the drag coefficient by Smith (1980).
(Fig. 8).This is the water temperature at which Henry's constant equals 1 for a 0 ppt salinity and 1 atm air pressure.The CO 2 concentrations in mol m −3 on the x-axis correspond to the CO 2 concentrations of 200 to 900 ppm in water at 17.38 • C and 0 ppt.Water temperature and salinity had a dual effect in the CO 2 flux across the air-water interface.The changes in the y-intercept were due to their effects in the solubility of CO 2 (k H ), whereas the steepness of the slopes was given by their effects in the water-side transfer velocities (Fig. 8a).The same test was done isolating the CO 2 term (Fig. 8b).Water temperature and salinity only affected the yintercept of the functions due to their effects in the solubility of CO 2 .All slopes exhibited the same steepness, as the transfer velocity was not included in the function.A fairly similar process occurred with the effects of air pressure.

Model application
The model was tested by comparing the CO 2 flux estimates with the CO 2 fluxes observed in Ria Formosa's main channels and at the nearby coastal ocean with the IRGA and floating chamber technique (Fig. 9).The model estimates were forced by the data on the environmental variables that were simultaneously collected.Data were not available to allow for estimates of roughness length inside Ria Formosa.Therefore, given the calm weather and smooth sea surface, these were arbitrarily given the value of z 0 = 10 −4 m (see Mackay and Yeun (1983) and Vickers and Mahrt (2006)).It tended to increase excessively the gas flux estimates on 3 March and 15 April.The stable and near neutral atmospheric conditions predicted for the April surveys did not have much effect on the gas flux predictions.On the other hand, the atmospherically unstable conditions predicted for 3 March overestimated the gas flux.For the nearby coastal ocean, the inclusion of sea-state and atmospheric stability was crucial for predictions to approximate the observations (in Fig. 9).Using the adapted Taylor and Yelland (2001) formulation with A = 1200, B = 4.5 and = 1 yielded roughness lengths around 10 −7 m to 10 −6 m and very poor fits (not shown).These improved significantly when = 0.355, A = 1.26 and B = 1.2 were used.Changing from the Mackay and Yeun (1983) (1975), 'Stu88 ' Stull (1988) and 'Le97' Lee (1997) based on 'B ' Businger et al. (1971) or 'D ' Dyer (1974).whitecap was fundamental at setting the water-side transfer velocity.The overall transfer velocity (i.e. the harmonic mean) estimated either from the water or air point of view (Fig. 10) was limited by the water-side transfer velocity.Their values are always approximate.On the other hand, the air-side transfer velocity was always about two orders of magnitude faster, proving it was never limiting the exchange.
The CO 2 in Ria Formosa's water body showed a pattern very similar to the CO 2 flux (Fig. 11) still, with a smoother variation.Here, positive values represent depletion (forcing uptake), whereas negative values represent surplus (forcing escape) of CO 2 in the water relative to what would be expected if it was in equilibrium with the overlying atmosphere.The heterogeneity was evident of the Ria Formosa water body in terms of CO 2 budget.In March it was behaving autotrophically, with a depletion of CO 2 relative to the atmosphere, whereas in April it showed an erratic behaviour, changing from autotrophic to heterotrophic in just a few hours.

Tuning the decomposition of the difference in the gas fluxes
The DDF analysis must be optimized before its application with the intention to minimize both the error in the estimates and its computational effort.This includes choosing for each of the tested variables (x i ) the order of the partial derivatives ( i ), the size (δ i ) and number (n i ) of the steps taken,  Mackay and Yeun (1983) and atmospheric stability scheme by Stull (1988).Wilke and Chang (1955), CD by Smith (1980) and k a by Mackay and Yeun (1983).and the point of estimation of the partial derivatives (k i ) in units of steps from x i,a .Knowing the computational effort and the error in the estimates are inversely proportional, it was searched for the right balance.The inference of the best options was summarized in Table 1 and Fig. 12 and 13.The cpu time was estimated for the n i steps in the tested variable with n i = 1 for all other variables.Not all possible model variables were tested but only the ones currently used for the  2001) and (ψ) accounting for atmospheric stability following Stull (1988) Wilke and Chang (1955), and k a by Mackay and Yeun (1983).
CO 2 flux estimates.The water temperature was set aside in Fig. 12 to exhibit a graphical representation of the typical evolution of the error.For this variable, as well as others like air temperature, salinity and wind speed, the optimal choices depended on the algorithms used.This work provides many optional algorithms, and it was not feasible to test them all.Only a few were tested and presented in the results.This does not mean these few were the best at estimating the CO 2 flux and should always be preferred.The optimization process also diverged whether the CO 2 concentrations were given in units of ppm or mol m −3 .Generally, using the mol m −3 units gave more accurate or equally accurate results and with less effort, the exception being with air temperature where it was the other way around.Fitting the gas flux to the z 0 , w and z using a multivariate collocation polynomial was only accurate if the n i steps of size δ i closely covered the range ( x i ) between the reference (x i,a ) and alternative (x i,b ) situations.It was attempted to feed the z 0 , w and z to the DDF tool in several units while adjusting the size and number of steps taken so that δ i would always be close to 1.There was no globally better solution.In the examples shown in Table 1, the best options were to give w in hm h −1 and taking 5 steps of size 3.6, z in dam (10 m) and taking 5 steps of size 1.29, and z 0 in mm and taking 1 step of size 8.9.
It was tested for the optimal point of estimation of the partial derivatives -that is, how far away from x i,a the partial derivatives could be estimated (Fig. 13).This distance is k i in the collocation polynomial and is given in units of steps taken away from x i,a .The k i need not be an integer number, as it was proved by testing it from 0 to 5 at 0.2 increments.The results are presented for the easy variable u 10 and three harsh variables T w , w and z.The partial derivatives of the easy variable could be equally well estimated at any k i within the bounds of x i,a and x i,b .On the contrary, the partial derivatives of the harsh variables could only be well estimated at k i = 0.
It was tested whether it is possible to estimate the collocation polynomial once for a hypervolume comprising all the desired samples and then estimate each sample-specific set of partial derivatives required for each DDF at its specific location within the hypervolume.This would relieve the software from estimating a new hypervolume for each new DDF, which is very time-consuming.On the available data set, it was considered x i,a as the minimum x i and x i,b as the maximum x i , over all samples.Then, the partial derivatives were estimated at point x i,c so that min x i ≤ x i,c ≤ max x i ,  inputting k i in units of steps of size δ i taken from min x i (i.e k i = (x i,c − min x i )/δ i ), and as long as δ i was always customized so that δ i • n i = max x i − min x i .The n 10 = 5 was important for the accuracy of the estimates of the partial derivatives related to w.The accuracy was generally remarkable (Fig. 14).Nevertheless, for the harsh variable of current velocity there were still a few cases for which they were very poor.This error was not due to the method being tested but rather due to the independent estimation of the partial derivatives, used for the comparison: whenever x i,c was to close to min x i or max x i , it forced δ i to be much smaller than 1, bringing severe error to these estimates.

Applying the decomposition of the difference in the gas fluxes
The decomposition of the difference between the CO 2 fluxes in the air-water interface inside Ria Formosa on 15 April 2011 for the first sample in the time series and in the nearby coastal ocean on 3 March 2011 (Fig. 15) had only a 0.04 % error relative to the CO 2 fluxes predicted by the model.This is the remainder of the Taylor expansion, i.e. the error specific to the DDF method.Still, it is known the flux predicted for 3 March was underestimated by about 2.5 mmol m −2 d −1 .Therefore, at least for a few variables their correct terms should be larger than the ones presented.The CO 2 flux was positive in the coastal ocean, meaning CO 2 uptake, whereas it was negative inside Ria Formosa, meaning CO 2 escape from the water to the atmosphere.As the coastal ocean was the reference situation (f a ), and Ria Formosa the alternative situation (f b ) the difference (f b − f a ) was negative.The biggest contributor to this difference was the CO 2 concentration in the water (C w ) as the coastal ocean was behaving autotrophically on 3 March and Ria Formosa was behaving heterotrophically on 15 April, at least at that section and between 11:00 and 13:00 LT.A smaller CO 2 concentration in the air (C a ) over Ria Formosa also gave a significant contribution to the CO 2 flux difference.The effect of atmospheric stability/instability was estimated aside from those of air (T a ) and water (T w ) temperatures.Then, the remaining effects of temperatures, salinity (S) and air pressure (P ) had only slight influence on the flux difference.The wind velocity (u 10 ) had a small negative term, because it was slightly windier on 15 April (4.5 m s −1 ) than on 3 March (3 m s −1 ).Nevertheless, the coastal ocean surface was much rougher (z 0 ) than the water surface at the lagoon system, generating more drag, a much higher k wind w and thus the large z 0 positive term.Also, the general wind transferred more momentum (relatively) to the air in contact with the coastal ocean surface given the atmospherically unstable conditions verified on 3 March 2011 than to the air in contact with the Ria Formosa water surface given the atmospherically stable conditions verified on 15 April 2011 -thus, the positive term for ψ u .It was assumed the misfit between the observed and predicted CO 2 fluxes at the coastal ocean was due to the underestimation of z 0 (due to uncertainty in the parameters) and/or ψ u (by neglecting cool skin and warm layer effects).Therefore, it was expected that the correct DDF terms for these variables to be larger.These positive terms discount from the overall negative sum, meaning that if the Ria Formosa on 15 April had its water surface as rough and the overlying atmosphere as unstable as the coastal ocean had on 3 March, the CO 2 flux difference would be even higher as there would be more transfer velocity and thus more CO 2 being transferred to the air over Ria Formosa.Nevertheless, inside the lagoon system the turbulence from below, that is from current drag with the www.ocean-sci.net/9/355/2013/Ocean Sci., 9, 355-375, 2013 bottom, compensated for the lesser turbulence from above, as it is shown by the negative terms relative to current velocity (w) and depth (z).

Model implementation
The application of the present model to estimate the CO 2 flux across the air-water interface showed the overall transfer velocity to be limited by the water-side transfer velocity.This is the expected for sparingly soluble gases such as CO 2 (Upstill-Goddard, 2006;Johnson, 2010).In this case the inclusion or not of the air-side transfer velocity and the choice of its formulation were irrelevant.The fundamental aspect was the water-side transfer velocity and the algorithms chosen to simulate it.On the contrary, for gases that are very soluble or react with water, the air-side transfer velocity is expected to be the limiting factor (Upstill-Goddard, 2006).In these cases the inclusion of the air-side transfer velocity should be crucial to accurately simulate the gas fluxes.Sander (1999) provides an extensive list of gases and their solubility in water.The estimation of the overall transfer velocity by the harmonic mean of the air-side and water-side transfer velocities weighted by the gas solubility proved to be an effective way to simulate this dynamics.Many different algorithms are available in the literature to estimate the water-side transfer velocities.The simpler ones are empirical formulations relating to the effect of a single factor such as wind, whitecap or current.Allowing for a variable drag coefficient dependent on wind speed, seasurface agitation and other physical properties of the atmo-  and Chang (1955), k a by Mackay and Yeun (1983) and ψ u by Stull (1988).(C a ) CO 2 in the air, (T a ) temperature, (P ) air pressure, (u 10 ) wind speed, (z 0 ) roughness length, (ψ u ) atmospheric stability (C CO 2 in the water, (T w ) water temperature, (S) salinity, (w) water current, (z) depth.Vectors n and were [2 2 1 2 3 3 1 3 1 2 2 (1958), α by Johnson (2010), D w by Wilke and Chang (1955), k a by Mackay and Yeun (1983) and ψ u by Stull (1988).(C a ) CO 2 in the air, (T a ) air temperature, (P ) air pressure, (u 10 ) wind speed, (z 0 ) roughness length, (ψ u ) atmospheric stability (C w ) CO 2 in the water, (T w ) water temperature, (S) salinity, (w) water current, (z) depth.Vectors n and were [2 2 1 2 3 3 1 3 1 2 2].
sphere and ocean (Sethuraman and Raynor, 1975;Smith, 1980;Mackay and Yeun, 1983;Smith et al., 1992;Taylor and Yelland, 2001) increases substantially the model accuracy.It is equally important to consider the advective components of k a and k w tend asymptotically to zero as the atmosphere changes to still air and the sea changes to still or deep water, the diffusive transport becoming the dominant feature.Therefore, any model parameterization meant to be applied to coastal studies and inland waters, where low wind is frequent, should account for the diffusive component of k a and k w and hence force them to stabilize in accurate values as turbulence decreases.However, most of the available formulations either neglect the diffusive transport or show great discordance about their related transfer velocities, revealing the lack of care this subject has been devoted.
Wind-based algorithms developed from open ocean data are usually second or higher order polynomials that increase the transfer velocity enormously with wind speed.Still, there is great variability within this set of algorithms.The wind based algorithms developed for coastal systems by Carini et al. (1996) and Borges et al. (2004b) are linear functions that underestimate the transfer velocities at high wind speeds relative to the open ocean formulations.Regarding this matter, three points should be taken into consideration: (i) at higher wind speeds the open ocean formulations are still interpolating, whereas the coastal system formulations are extrapolating; (ii) at wind speeds as high as 30 m s −1 , even the open ocean formulations are extrapolating; and (iii) fetch is a fundamental aspect not taken into consideration in any of these formulations and exhibits its widest change precisely when compared between coastal systems and open ocean.Raymond and Cole (2001) fit an exponential function to data from estuaries collected at low winds.Extrapolation to high winds yielded transfer velocities outstandingly higher than any other, even for open ocean.This is probably the best demonstration that the application of many transfer velocity algorithms should be restricted to the specified environmental conditions upon which they were developed.
Slightly more elaborated algorithms integrate the effects of a few factors, allowing for an increase in their applicability and accuracy.However, most of these are still empirical relations constrained to the environmental range upon which they were tested.Considering the broad applicability to the coastal ocean, rivers, estuaries and lagoonary systems, it is relevant that only the numerical schemes by Borges et al. (2004b) and Johnson (2010) comprise the effect of salinity changes and only the one by Borges et al. (2004b) imports the effects of current drag from previous authors.A few numerical schemes have gone further with more mechanistic approaches to the environmental processes they represent.This allows for a significant increase in their applicable environmental range and possible interaction with complementary formulations.It is the particular case of Memery and Merlivat (1985), and Johnson (2010) that the COARE algorithm and the vast body of literature are related to the surface renewal theory and micro-scale wave breaking.
The present numerical scheme tries to incorporate all these options and develop a software able to estimate the gas flux across the air-water interface under the broadest range of environmental conditions with a unique model parameterization.The estimates of the transfer velocity showed that in shallow coastal waters the effect of water current can be as important as the effect of low to moderate winds.In macro-and mesotidal estuarine and lagoonary systems, higher tidally driven water currents occur on a daily basis, whereas high winds do not.Therefore, the effects of water current and depth are fundamental for the model performance in coastal environments.On the other hand, the attempts to calibrate the model for the coastal ocean samples demonstrated the importance of roughness length and atmospheric stability for the estimation of the gas fluxes across the surface of large water bodies.
The roughness length formulation by Taylor and Yelland (2001) is very practical as it requires only two parameters from the wave field.It is also very intuitive as it states the roughness length scaled to the wave height is proportional to the wave slope, this function being linear or exponential depending on the exponent (B) value.However, the wave fields are not uniform and may be decomposed into a wave spectrum where each of its components potentially gives a relative contribution to the roughness length.The alternative proposed by Taylor and Yelland (2001) is to use the peak component of the wave spectrum.This simplification may imply loss of information and predictive power.However, Moon et al. (2004) have demonstrated for tropical cyclones the Charnock coefficient is mainly determined by wind speed and the peak wave age, thus supporting such simplification.This problem is aggravated by the fact that roughness length is a theoretical concept that cannot be tested directly.Usually are used proxies such as the friction velocity, the drag coefficient or the Reynolds number (Sethuraman and Raynor, 1975;Taylor and Yelland, 2001;Fairall et al., 2003;Frew et al., 2004;Moon et al., 2004).
The present simulations have demonstrated atmospheric stability to have a huge potential to influence the friction velocity and therefore the transfer velocity.Nevertheless, the evaluation of atmospheric stability and its application to marine coastal environments should be cautious as Vickers and Mahrt (2006) propose that (i) Monin-Obukhov similarity theory does not apply to sea surfaces with sharp temperature gradients and (ii) the sensible heat flux is better correlated with the sea surface temperature in a 1-2 km downstream lag.Here may lie the explanation for the predicted atmospherically unstable conditions overestimating the CO 2 flux on the 3 of March survey inside Ria Formosa.At the coastal ocean the atmosphere boundary layer was only in contact with the sea surface, and thus their temperature gradient was estimating real atmospherically unstable conditions.Inside Ria Formosa was not high tide, and the atmosphere boundary layer was in contact with many other surfaces besides the sea surface.Their temperature gradient was estimating probably unreal atmospheric unstable conditions.The formulation by Stull (1988) is apparently the most reliable.The formulation by Lee (1997) had two problems.The first was it did not work for stable conditions of Ri b > 0.2.This is a problem common to many other formulations based on the Businger et al. (1971) and Dyer (1974) works.The other was that for very unstable conditions (Ri b ≤ −0.2) the function started behaving exponentially, which is unreal.In the publication by Lee (1997) this does not happen.Yet, the algorithm was scrupulously imported to the present model and software.Both these formulations rely on the famous works by Businger et al. (1971) and Dyer (1974).These had robust experimental designs and data sets, becoming the core of the science about atmospheric stability.The work by Sethuraman and Raynor (1975) was more modest.The small data set did not have | T |> 3 • C and was very scattered around the linear fit, bringing legitimate suspicions about the adequacy of the equation structure and parameter values.When T > 3.7 • C, the equation predicts a negative drag coefficient, which is physically impossible.When T < −3 • C, the equation quickly tends to predict hurricane force friction velocities.
Under high winds the effects of whitecap and bubbles become important (Memery and Merlivat, 1985;Zhao et al., 2003;Woolf, 2005) and therefore should be added to the model.Memery and Merlivat (1985) propose a complex algorithm that accounts for many physical properties of water and bubbles.Woolf (2005) states the water-side transfer velocity as the addition of a term for the breaking waves and another for non-breaking waves.Presently it was only implemented the simpler solution by Zhao et al. (2003).It is debatable whether this formulation should be overlapped or not with the roughness length formulation as ultimately both account for the effect of the wave field in k wind w and therefore may be redundant.In the preliminary test performed in this work, their overlap gave the best results, far beyond any other.
In the coastal ocean, as the swell approaches the shore, the drag with the shallower bottom compacts the waves, decreasing the wave length while keeping the wave height.The wave slope increases and thus also the roughness length (Taylor and Yelland, 2001).Therefore, these authors expect the gas transfer velocity to increase as the coastal ocean approaches the shore.This has two implications for the current work.One is that the data from field surveys or oceanographic numerical laboratories should not neglect the effect of increasing wave slope with decreasing depth.The other is to clarify that the surface roughness in Ria Formosa is generated exclusively inside the lagoon and independent from the swell outside.Nevertheless, it should be considered the possibility that the downwind depth profile inside estuarine and lagoonary systems may have an effect in roughness length and consequently in the gas transfer velocity, as Upstill-Goddard (2006) proposes for generalized shallow waters.Also the presence of surfactants decreases the gas transfer velocity (Memery and Merlivat, 1985;Frew et al., 2004), particularly with lower wind speeds, and surfaces with shorter waves are more affected by surfactants (Frew et al., 2004).Therefore, a likelier presence of surfactants inside estuaries and lagoons than in the nearby coastal oceans should also be considered.
Finally, the current software allows for the gas concentrations to be input in units of ppm, although the model requires them to be converted to units of mol m −3 .This conversion is dependent on temperature, pressure and salinity, and thus is yet another way to account for the effects of these variables in the flux of a gas across the air-water interface.This is not a model artificialization but rather represents simple objective environmental features.Taking the example of the atmosphere, as an air mass changes its density, it keeps its inner relative gas concentrations (given in ppm) but changes its volumetric gas concentrations (given in mol m −3 ), thus affecting its gas exchanges with any other distinct entity.

Model alternatives
The quantification of the effects of wind and surface roughness was done following two alternatives: (i) the wind loglinear profile and (ii) the Sethuraman and Raynor (1975) formulation.It is also possible to use the Frew et al. (2004) transfer velocity formulation where it is the exponent upon the Schmidt number to show a dependency on sea surface roughness.For an indirect estimation of roughness, length from the wave field was used Taylor and Yelland (2001) formulation relating surface roughness to the wave slope.Frew et al. (2004) and Hwang (2005) present alternative formulations based on the mean square slope.The wave slope may be estimated using a pressure transducer.However, in the smooth surfaces that often occur in estuarine and lagoonary systems under calm weather are required pressure transducers with resolutions higher than 4 Hz.Alternatively, the wave field may be estimated using a scanning laser slope gauge (Frew et al., 2004).
Alternatively, friction velocity may be estimated from the roughness length following Charnock's model (Charnock, 1955): z 0 = α c u 2 * /g, where α c is Charnock's coefficient.This alternative implies estimating α c as a function of u 10 and the input wave age (Moon et al., 2004).Friction velocity may also be estimated from the near surface covariance of horizontal (u ) and vertical (w ) wind components.Then, the gas flux model and DDF analysis must account for friction velocity directly and in replacement of roughness length (z 0 ) and wind speed (u 10 ).For the model estimation the calculus is simpler as it is a simple function of the horizontal and vertical variability of the wind components.Nevertheless, as for all the alternatives presented that require simpler calculus, these have the cost of information being lost for the DDF analysis.For the example shown in this work, it would not be possible to assess whether (or how much of) the difference between the CO 2 flux inside the lagoon and in the coastal ocean was due to the difference in the wind properties or due to the difference in the sea surface roughness.
The total transfer velocity of a gas may also be estimated from the total transfer velocity of heat (Frew et al., 2004).The relation is given by k gas = k heat (Sc/P r) −n , where Sc is the Schmidt number, P r the Prandtl number and n a scalar (usually between 0.5 and 0.7).In its turn k heat = j heat /(ρc p T ), where ρ and c p are seawater density and specific heat, respectively; T is the seawater temperature difference between the "cool skin" and the bulk of the surface boundary layer, which may be estimated from infrared imagery; and j heat is the net heat flux density at the sea surface, which may be estimated from micrometeorological measurements.

Tuning the decomposition of the difference in the gas fluxes
When performing the DDF, it is intended to have the most accurate results.Still, not wasting time, generally it is not worth taking more steps than the optimal order of the partial derivative.The optimal choices varied with the numerical options but also with the units used to give the CO 2 concentrations.This latter was because several environmental variables affected the solubility/volatility and therefore the conversion of the CO 2 concentrations when given in ppm to the mol m −3 units required by the flux model.
The DDF optimization relative to z 0 (x 5 ), ψ u (x 6 ), w (x 10 ) and z (x 11 ) was more complicated, because fitting the gas flux using an n-th order collocation polynomial was only accurate if the n i steps closely covered the h i range (for i equal to 5, 6, 10 and 11).This obliges conjugating h i with (1) the chosen x i units to feed the model, (2) the n i steps taken and (3) the δ i size of the steps taken.Therefore, the optimization of the DDF relative to these variables must always be customized to the data set.A good rule of thumb is to choose the units so that h i has one digit.Afterwards, δ i should equal the maximum h i found for all alternative sites divided by n i keeping in mind that δ i should never become too big nor too small.In order to illustrate the relevance of such a procedure, the optimization of the DDF relative to the depth parameterization was intentionally shown (in Table 1) for a reference site at open ocean (z = 67 m) and an alternative site inside Ria Formosa (z = 2.5 m).Big depth differences may occur in future applications of this DDF tool.Therefore, it was essential to show that the collocation polynomial is so sensitive to depth that n i •δ i must match h i for the DDF to be accurate.In this case it was 5 1.29 dam =6.45 dam.However, this DDF tool is also intended to be applied to several (possibly many) alternative situations, and it is not practical for the user to have to customize δ i by hand for each new alternative situation.Therefore, the software was updated to do it automatically, whenever required by the user, to whatever variables selected, in whatever units fed to the DDF tool, by setting δ i = h i /n i .With this customization may occur a hidden bias passing undetected.When n i •δ i is very close to h i , the Taylor series always closely matches f b −f a , irrespective of n i .This implies that the estimated error (1 minus the sum of all the terms) is very low although each term individually may be biased -in fact, even if relevant higher order terms are missing.The end result is a very low estimated error although the partition of f b − f a among the several environmental variables is severely biased.To overcome this problem the choice of i and n i must be independent of this customization process where δ i = h i /n i .
One important and immediate application of this DDF tool is to conjugate it with numerical modelling labs such as MO-HID, ECO lab, URI's, WRL's or FIO's.These numerical labs simulate the evolution of the physical, chemical and biological properties of the marine and aquatic environments in a particular area.In order to do that the domain area is often divided into thousands to tens of thousands of smaller units.The evolution of the model properties is often estimated at time intervals of a few seconds.It is unfeasible to apply the DDF tool to thousands of locations every few seconds.However, it is possible to drastically lighten it up to the point of enabling this application.The feature that turns the DDF algorithm computationally heavy is the estimation of the hypervolume of multivariate finite differences needed for the estimation of the partial derivatives.The execution of this calculus for each point and iteration is what makes its application unbearable.The solution to this problem relies on the estimation of this hypervolume only once for the whole spatial domain over a major time interval.This hypervolume must then comprise a grid that, for each of the environmental variables, stretches from the minimum to the maximum recorded values, including reference and all alternative sites.Afterwards, it is possible to accurately estimate the partial derivatives at any point inside this grid, because the algorithm used for its estimation (presented in Supplement B) works equally well for k i being an integer or fractional number.The tests to the estimation of the partial derivatives at any point x i,c inside this grid gave a remarkable accuracy, proving this to be the right solution.

Insights into the subject system
The gas flux models integrated with the DDF have shown to be valuable tools for the study of any gas crossing the airwater interface, may it be a pollutant or part of a biogeochemical cycle.The gas flux numerical scheme allows choosing the empirical formulations most suited to a particular case or, alternatively, mechanistical formulations of broader application.It further allows identifying past cases where inappropriate parameterizations may have been used and quantifying the expected biases.As an example, Oliveira et al. (2012) studied the Portuguese coast as a sink/source of CO 2 .For that they estimated its flux between the atmosphere and the coastal ocean adjacent to the Douro, Tagus and Sado estuaries.The fluxes were estimated from the formulations by Carini et al. (1996), Raymond and Cole (2001) and Borges et al. (2004b) applied to measures of the required environmental variables.However, actual field measurements of the fluxes were not done, which would enable validation.The problem here was that these formulations were neither developed from open ocean data nor supported by data on high wind conditions.While the use of the Carini et al. (1996) and Borges et al. (2004b) parameterizations clearly underestimates the flux at open ocean, the extrapolation of the exponential function by Raymond and Cole (2001) is a very wild guess.To illustrate it, during the cold front rampant over Europe, the water off-shore Ria Formosa on 4 February 2012 by 09:50 LT was at 15.1 • C, the significant wave height at 1.54 m, the wave length at 31.6 m, the average wave period 4.5 s, the air at 6 • C and the wind blowing off-shore at 10 m s −1 .Following the Taylor and Yelland (2001) formulation with its original parameters, z 0 =2.3 mm.Following Stull (1988), Ri b = −0.04,ψ u = −0.81 and u * =0.53 m s −1 .Given these conditions the sea-surface roughness, whitecap and atmospheric instability should play a major role setting the water-side transfer velocity.These estimated by the Carini et al. (1996) and Borges et al. (2004b) formulations are of 20.3 and 26.8 cm h −1 , respectively.The estimated by Raymond et al. (2000) formulation is of 63.3 cm h −1 .When estimated from wind log-linear profile and the Zhao et al. (2003) formulation for the effect of wind and whitecap, it is of 56.2 cm h −1 .But, if atmospheric stability is estimated from the formulation by Lee (1997) based on Businger et al. (1971), the u * = 0.64 m s −1 and the water-side transfer velocity is 71.6 cm h −1 .When using the CO 2 flux across the air-water interface as a proxy for the ecosystem metabolism, one must take into account that it is also strongly dependent on the influence of turbulence on the transfer velocity.To correct for this, Frankignoulle (1988), Smith and Hollibaugh (1993), Raymond et al. (2000), Cole and Caraco (2001), Koné et al. (2009) and Torres et al. (2011) tested using only the difference term of the flux equation.By decomposing the fluxes in all their parcels, the DDF further allows accurate estimates of the influence of a wide range of environmental variables in mediating the flux, together with its spatial and temporal variability.Furthermore, this tool allows focusing on the effect of a specific variable at different places, different times or under different methodologies filtering out the undesired effects of changes in other variables.It is also possible to use the DDF tool focusing on a specific aspect.For the subject of the drives for a transfer velocity, it only requires replacing the flux by the transfer velocity as the dependent variable.Similarly, for the subject of the drives for a difference between the gas concentrations in the air and water phases, it only requires the replacement of the flux by the C or ppm as the dependent variables.
The gains brought by this new gas flux numerical model and DDF tool were clearly demonstrated with the comparison between Ria Formosa on 14 April 2011 and its surrounding coastal ocean on 3 March 2011.While the coastal ocean was behaving autotrophically, the Ria Formosa was behaving heterotrophically, at least between 11:00 and 13:00 LT and at that particular site.The bulk of the CO 2 flux difference was indeed due to the difference in the CO 2 concentrations in the water inside and outside.However, there were also other factors taking part that the DDF enabled to set aside.While the transfer velocity in the ocean was set by turbulence from above, inside the mesotidal lagoon system it was majorly set by turbulence from below.A similar contrast was presented by Borges et al. (2004a) when comparing between micro-, meso-and macrotidal estuaries.On the other hand, Ho et al. (2011) determined the transfer velocity in the Hudson River was basically set by wind speed and independent of current drag with the bottom.Still, these authors admitted such results may have been influenced by samples having been taken tendentiously over ≥ 5 m depths.
The CO 2 series suggest Ria Formosa could have been behaving autotrophically in early March when the water was around 17 • C and could have been behaving heterotrophically in mid-April when the water was around 20.5 • C.This change with temperature may be related to the dominant biological process taking place.Photosynthesis by seagrass meadows is much less sensitive to temperature changes than respiration by bacteria.Furthermore, the water column at the sampled Ria Formosa channel during ebb tide changed from autotrophic to heterotrophic in a couple of hours, proving there was a strong spatial/temporal heterogeneity in the CO 2 balance.It is hypothesized the metabolic status of a particular section of the water column was related to it being over a seagrass meadow or a mudflat in its near past.Seasonal and short-term shifts of the CO 2 balance in estuaries, lagoonary systems and coastal waters were already reported by Raymond et al. (2000), Cole and Caraco (2001), Frankignoulle et al. (2001), Borges et al. (2004a), Borges (2005), Koné et al. (2009), Hunt et al. (2011), Torres et al. (2011) and Oliveira et al. (2012).

Conclusions
Wide spatial and temporal variabilities of gas concentrations in the water, in the overlying air and their fluxes across the air-water interface are widely documented for the open oceans, the coastal oceans and riverine systems.These gas fluxes have a multitude of potential forcing functions.However, their integration and the establishment of their relative importance has been underachieved.This is particularly evident from how atmospheric stability, sea-surface roughness and current drag with the bottom have often been devalued in studies about riverine systems and coastal waters.The currently presented numerical tools give a significant contribution to this subject.Now it is easier to use a single model for any type of marine and freshwater environment and to conclude the differences found between those report exclusively to the environments and not to different numerical options.Furthermore, the numerical scheme allows for the upgrade of each relevant environmental process already implemented as well as the addition of new processes.Any interested researcher is free to add a particular formulation for his/her own personal use and is further invited to share it with everyone else.The versatility of the present model, tools and software allows the user to follow two distinct approaches.The user may choose to use the formulations available in the literature that best fit to a particular situation.These tend to be more of an empirical nature and to fail under largely different environmental conditions.Alternatively, the user may build the model upon a more mechanistic approach, computationally heavier, but tending to yield better global fits.The DDF tool allows for the quantification of the effects of all the environmental variables and processes involved in the gas flux across a particular air-water interface relative to a reference one.It further allows focusing on a specific variable or process eliminating the error from the remaining ones.

Fig. 7 .Fig. 7 .
Fig. 7. Effect of surface roughness (measured by z 0 ) and atmospheric stability (measured by ∆T = T a − T w ) in (a) the air-side and (b) the water-side transfer velocities.T w = 15 • C, u 10 = 5 m s −1 , k a and k wind w byMackay and Yeun (1983) and atmospheric stability scheme byStull (1988).

Fig. 10 .Fig. 10 .
Fig. 10.Transfer velocity limiting phase.Overall transfer velocity from the air (k(a)) and from the water point of view (k(w)); air-side transfer velocity (k a ) and water-side transfer velocity (k w ).

Fig. 12 .
Fig. 12. Choosing n and Θ for the T • w C terms in the DDF.(a) n = Θ and (b) n = 5.The (a) and (b) estimates with [CO 2 ] given in ppm.(c) n = Θ with [CO 2 ] given in mol m −3 .

Fig. 12 .
Fig. 12. Choosing n and for the T • w C terms in the DDF.n = and (b) n = 5.The (a) and b) estimates with [CO 2 ] given in ppm.(c) n = with [CO 2 ] given in mol m −3 .

Fig. 14 .Fig. 14 .
Fig.14.Accuracy of the partial derivatives estimated at k i located between min x i and max x i .Results are shown for wind speed (u 10 ), water temperature (T w ) and current velocity (w).Derivatives estimated by (f) forward formula and (b) backward formula.

Fig. 15 .
Fig. 15.Decomposition of the difference between the CO 2 Flux in Ria Formosa on the 15 April 20 and in the nearby coastal ocean on the 3 March 2011.Transfer velocity (k) by double layer, k wind w

]. 55 Fig. 15 .
Fig. 15.Decomposition of the between the CO 2 Flux in Ria Formosa on 15 April 2011 and in the nearby coastal ocean on 3 March 2011.Transfer velocity (k) by double layer, k wind w

Table 1 .
Optimization of the DDF.* any D w scheme.(ad.fit.) δ i adjusted to fit h i .Optimality is bolded.