Gas Transfer Velocity from QuikSCAT/SeaWinds Backscatter

A major goal of NASA's Earth Science Enterprise is to assess the net global flux of CO₂, a climatologically important, radiatively active gas, into the ocean (Asrar and Dozier, 1994). Flux estimates based on air-sea exchange models are an important diagnostic tool and must be consistent with the long-term sequestering of CO₂ predicted by various carbon cycle models, both general circulation models (GCM) and isotopic estimates (¹⁴C). In certain regions of the ocean, exchange across the air-sea interface may be the rate-limiting step, as in high latitudes where episodic deep advection processes occur on time scales that are shorter than ocean-atmosphere equilibration times (Broecker and Peng, 1982). These regions apparently account for much of the sensitivity of carbon flux models to variations in the gas transfer coefficient (Sarmiento et al., 1992; Johnson, 1995).

Gas exchange flux is calculated from Eqn. (1) which contains two factors: a concentration gradient across the air-water interface and the exchange coefficient or transfer velocity representing a parameterization of a combination of important near surface exchange mechanisms.

(1)

Here F is the flux of mass (mol cm^-2 hr^-1), k is the transfer velocity (cm hr^-1) and Δ C is the concentration difference (mol cm^-3). Transfer velocity fields predicted by various parameterizations based on wind speed (Liss and Merlivat, 1986; Wanninkhof, 1992; Tans et al., 1990; Erikson, 1993) lead to widely varying estimates of zonal and global net CO₂ fluxes (Etcheto and Merlivat, 1988; Boutin and Etcheto, 1995). These are not sufficiently constrained to validate GCM models (Sarmiento et al., 1992; Keeling et al., 1989; Stocker et al., 1994), which suggest a global uptake of 2±0.8 GtC/yr (in the mid-1980's), or to shed light on the apparent ``missing sink'' for anthropogenic CO₂ (~1.6 GtC/yr). The uncertainty is a significant fraction of the total annual 3.5 GtC uptake by non-atmospheric sinks (Johnson, 1995). Thus, the Intergovernmental Panel on Climate Change (IPCC, 1996) has identified uncertainty in the gas exchange coefficient as a significant limitation in assessing the role of the ocean in absorbing anthropogenic CO₂ and has called for increased study of its global spatial and temporal variations in order to help close the global carbon budget.

Current parameterizations of k are based on wind speed at 10 m height, U₁₀ (Liss and Merlivat, 1986; Wanninkhof, 1992; Wanninkhof and McGillis, 1999; Nightingale et al., 2000). Such parameterizations would be extremely useful in estimating CO₂ fluxes both seasonally and spatially, since global wind fields can be estimated from space-based scatterometers (Chelton et al., 1990). Recent studies, however, have shown that k is not a unique function of wind speed (Frew et al., 1995; Frew, 1997; Hara et al., 1995, Bock et al., 1999). Evidence for this is seen in the considerable scatter in both laboratory and field data when correlated with either U₁₀ or friction velocity u_* (Wanninkhof, 1992). Other factors have a strong influence on k, most notably wave fetch (Wanninkhof, 1992), boundary layer stability (Erickson, 1993), and the presence of surface-active organic matter, which affects the small-scale wave field and surface turbulence (Frew et al., 1990; 1995; Bock et al., 1995; Frew, 1997). These factors are expressed largely as modulations of surface roughness (and hence k).

We have avoided complications introduced by the U₁₀ model by developing an alternative method for predicting k using the TOPEX dual-frequency normalized altimeter backscatter (Frew et al., 2005; Glover et al., 2002). The modulating factors cited above are assimilated using a direct measure of surface roughness, the mean square surface slope, which has been shown to be directly related to the transfer velocity (Jähne et al., 1987; Hara et al., 1995; Bock et al., 1999). The transfer velocity is derived from the inverse relationship between σ^o and the mean square slope () of the wave field from which it is reflected (Jackson et al., 1992). Winds derived from altimeters, radiometers and scatterometers are all empirical estimates based on surface roughness scaled to obtain wind speed. Nonetheless, the relation between U₁₀ and k cannot be specified precisely. The use of altimeter data instead of other sensors such as SSM/I or the NSCAT replacement (SeaWinds), which would provide better spatio-temporal coverage was attractive because of the existence of a rationale (theoretical model and empirical data) that allowed us to relate roughness to gas exchange rate. This report details the status of our project to extend our altimeter-based algorithm to a scatterometer-based algorithm by using altimeter results to bootstrap calibrate the scatterometer k-σ^o function.

We believe a radar backscatter-gas transfer velocity relationship represents a significant improvement over using wind speed to predict transfer velocity. While the initial development of the gas exchange algorithm has taken place using altimetry data (as members of the Jason-1 SWT), a relationship between scatterometer backscatter and transfer velocity has greatly accelerate further development of both algorithms and ultimately leads to a better data product in terms of both spatial and temporal coverage. Our contribution brings together scatterometer and altimeter remote sensing with in situ data from field programs. The combination of transfer velocity fields with pCO₂ fields derived from general circulation models (Gent et al., 1998) or field programs (Takahashi et al., 1997) allows us to identify the major CO₂ source and sink regions in the world's oceans and provides an unprecedented continuous record (using the combination of TOPEX/Poseidon, extended TOPEX/Poseidon, Jason-1, ALT, QuikSCAT and ADEOS-2 SeaWinds missions) of estimated CO₂ flux on monthly, annual, and decadal time scales.

2.0 Methods

Additional field measurements used to optimize these parameters were U₁₀ and a differential colored organic matter measurement (δCDOM). The field measurements of U₁₀ were made with a meteorological package onboard a research catamaran, adjusted to a neutrally stable 10 m height above the air-sea interface. The δCDOM were made with a fluorometer and represent the excess CDOM at the air-sea interface with respect to the bulk water concentration. These δCDOM measurements acted as a proxy for the presence of surfactants (Frew et al., 1990).

2.2 Match-ups

For each day, on a swath-by-swath basis, we map out all of the TOPEX 1 Hz radar returns that fall inside a QuikSCAT wind vector cell (WVC) restricted to be to within ±30 min of each other. A nonlinear least squares fit is performed with the TOPEX σ^o providing a

estimate on the left hand side and QuikSCAT providing the σ^o and φ information of the right hand side of Eqn 13. The parameters (p) are then used to process all QuikSCAT σ^o into an estimate of

, which is then used with the quadratic function that relates

to k (Eqn 20). The transfer velocities are computed at each wind vector cell (WVC) location in each orbital swath. Each of the level 2B WVC files are matched with the level 2A geolocated σ^o files. The selected wind speed and direction (DIRTH nudging algorithm) are extracted with the WVC quality flag (JPL, 2000). Only the wind vectors with a quality flag of zero are retained for further processing. All of the σ^o values for each valid wind vector are extracted along with the azimuth and incidence angles along with the 2^nd order polynomial coefficients for the weighting factor. The transfer velocity is then computed for each weighted σ^o within a WVC and averaged. The resulting averaged k is then saved with the latitude and longitude of the WVC.

2.3 Algorithm

The basis function for our algorithm is derived from two basic equations taken from the literature. From Apel (1994) we have:

(2)

where

is the mean square slope of the small scale wind-waves; U₁₀ is the wind speed at 10 m height; and α is a ``well known'' constant (Apel, 1994) equal to 0.004. In addition, we use the generic scatterometer model function from Stewart (1985):

(3)

where σ^o is the normalized radar backscatter, in this case when specifically from QuikSCAT written as; n(θ) an exponent that is a function of the incidence angle (θ) only; φ the relative azimuth, that is the difference between the direction the scatterometer is pointing and the direction the wind is blowing; and a, b, and c are constant parameters. While Eqn 3 is generally thought of as "the" scatterometer model function (Stewart, 1985), we note that Freilich points out in the QuikSCAT Algorithm Theoretical Basis Document (ATBD) that in reality the model function is much more complicated and therefore QuikSCAT does not use this simple function as a basis (http://eosdatainfo.gsfc.nasa.gov/eosdata/quikscat/quikscat_dataprod. html).

Since QuikSCAT has only two possible incidence angles (θ=46^o or 54^o), we can remove the dependence on θ by separating the data by incidence angle into two groups before proceeding with any further calculation. In this manner Eqn 3 becomes:

(4)

where n can now be treated as a constant parameter. Rearranging Eqn 2 as:

(5)

and directly substituting yields;

(6)

Now we can combine the constant parameters, such as:

(7)

and

(8)

Now we write the basis function for relating the mean square slope of the small-scale waves to the normalized radar backscatter from a scatterometer as:

(9)

where we represent σ^o in ``natural units'' (Freilich and Vanhoff, 2003), not dB i.e.,

(10)

This transformation is warranted and required by the non-linear least squares fitting routine being used to estimate q_{1, ..., 4} (Levenberg-Marquardt method, Press et al. 1992).

There is one additional, practical consideration to be made with respect to this algorithm: the non-negativity of . By definition is always positive and in the process of performing a non-linear regression this is easily assured by adding a small "DC-offset". The basis function used to fit the constant parameters to the matched TOPEX/QuikSCAT data would be:

(11)

From inspection, a non-linear regression using Eqn 11 as its basis function will result in parameters (q_{1, ..., 5}) composed of complex numbers. The reasons for this are: a) the factor is not always positive and b) the random noise found in the data prevents q₂ from always being an integer. Since it is not practical to use an algorithm whose operating parameters are complex numbers, this presents a problem.

The same inspection process does, however, suggest an alternate basis function that should both fit the data well and avoid complex number parameters. Note that Eqn 11 contains many similarities to Eqn 4. It appears that the function describing is a combination of a term that represents power-law relationship between and σ^o modulated by a double cosine function of the relative azimuth (φ). First we rename the parameters to p_{1, ...,
5} to make explicit this is a different basis function. Then by applying the second parameter (p₂) to only the σ^o variable (as ) and moving the part of the function into the numerator to act as a modulator, one arrives at the following basis function:

(12)

This function, however, is purely empirical in nature. Examination as to why this relationship works better than Eqn 11 is deferred to the Discussion Section.

There is still that one additional consideration. As in Freilich and Vanhoff (2003), Eqn 12 is valid only if σ^o is expressed in engineering or natural units. Since the QuikSCAT σ^o are report in units of dB, a transformation is required (Eqn 10) and Eqn 12 becomes:

(13)

Initially, this function was fit without the use of weighted . Equation 13 becomes:

(14)

where w is the measurement error-based weight of each obtained from the level-2A data second order polynomial expansion in σ^o, which expresses the variance (σ²) in the σ^o measurements, transformed to natural units.

(15)