Abstract

This study introduces exchange coefficients for wind stress (CD), latent heat flux (CL), and sensible heat flux (CS) over the global ocean. They are obtained from the state-of-the-art Coupled Ocean–Atmosphere Response Experiment (COARE) bulk algorithm (version 3.0). Using the exchange coefficients from this bulk scheme, CD, CL, and CS are then expressed as simple polynomial functions of air–sea temperature difference (TaTs)—where air temperature (Ta) is at 10 m, wind speed (Va) is at 10 m, and relative humidity (RH) is at the air–sea interface—to parameterize stability. The advantage of using polynomial-based exchange coefficients is that they do not require any iterations for stability. In addition, they agree with results from the COARE algorithm but at ≈5 times lower computation cost, an advantage that is particularly needed for ocean general circulation models (OGCMs) and climate models running at high horizontal resolution and short time steps. The effects of any water vapor flux in calculating the exchange coefficients are taken into account in the polynomial functions, a feature that is especially important at low wind speeds (e.g., Va < 5 m s−1) because air–sea mixing ratio difference can have a major effect on the stability, particularly in tropical regions. Analyses of exchange coefficients demonstrate the fact that water vapor can have substantial impact on air–sea exchange coefficients at low wind speeds. An example application of the exchange coefficients from the polynomial approach is the recalculation of climatological mean wind stress magnitude from 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40) data in the North Pacific Ocean over 1979–2002. Using ECMWF 10-m winds and the authors’ methodology provides accurate surface stresses while largely eliminating the orographically induced Gibb’s waves found in the original ERA-40 surface wind stresses. These can have a large amplitude near mountainous regions and can extend far into the ocean interior. This study introduces exchange coefficients of air–sea fluxes, which are applicable to the wide range of conditions occurring over the global ocean, including the air–sea stability differences across the Gulf Stream and Kuroshio, regions which have been the subject of many climate model studies. This versatility results because CD, CL, and CS are determined for Va values of 1 to 40 m s−1, (TaTs), intervals of −8° to 7°C, and RH values of 0% to 100%. Exchange coefficients presented here are called the Naval Research Laboratory (NRL) Air–Sea Exchange Coefficients (NASEC) and they are suitable for a wide range of air–sea interaction studies and model applications.

1. Introduction

Previously, Kara et al. (2000, 2002) introduced polynomial functions for determining the exchange coefficients of wind stress (CD), latent heat flux (CL), and sensible heat flux (CS) used in the bulk formulas. In these algorithms, the stability correction was approximated by a term depending on the sea–air temperature difference (TsTa), where Ta is at 10 m, and a second-order polynomial function of the wind speed at 10 m above the sea surface (Va) for CD or the inverse wind speed for CL. The algorithms do not require iteration to account for stability; therefore, they are computationally efficient for use in high-resolution coupled atmosphere–ocean circulation models and ocean general circulation models (OGCMs). They are also useful when, every model time step, latent and sensible heat are calculated via the bulk formulas using both ocean model sea surface temperature (SST) and the high-frequency (e.g., 6-hourly or daily) wind speeds. This approach keeps model SST on track without any explicit relaxation to SST data or a climatology (e.g., Wallcraft et al. 2003).

The algorithms presented in the two studies mentioned above were obtained using tabulated exchange coefficients, which were published by Bunker (1976) and Isemer et al. (1989). The use of these two datasets resulted in a few limitations in the algorithms. For example, the datasets did not provide realistic exchange coefficients at high and low wind speeds, predominantly to allow for the biases inherent in estimating fluxes from ship observations. Many later studies indicated lower values for the exchange coefficients (e.g., DeCosmo et al. 1996; Yelland et al. 1998, 2002), in general. As a result, the algorithms based on these two sets usually resulted in errors (5%–25%) in estimating the wind stress and latent and sensible heat fluxes in comparison to those obtained from the standard Coupled Ocean–Atmosphere Response Experiment (COARE) version-2.5b (v2.5b) algorithm (Fairall et al. 1996). The COARE v2.5b algorithm determines the turbulent air–sea fluxes through a recursive solution for the atmospheric dynamic stability. However, this algorithm also has some deficiencies because it was only verified for winds up to 12 m s−1. It is also noted that when Kara et al. (2000) was published and Kara et al. (2002) was under review, the COARE v2.5b algorithm was the accepted standard for air–sea flux estimates in the peer-reviewed scientific literature.

Fairall et al. (2003) later developed a new version of the COARE algorithm (v2.6 in 2001, renamed as v3.0 in 2003) with the aim of providing more realistic exchange coefficients at higher wind speeds based on the data collected since the release of COARE v2.5b. The new algorithm also takes additional observational data (Yelland et al. 1998; Hare et al. 1999) into account, which suggests that a Charnock constant value of 0.011 (Smith 1988), used in the earlier COARE v2.5b algorithm, is too low for higher wind speeds. Thus, a linear increase in the Charnock parameter is introduced between 10 and 18 m s−1 and a value of 0.018 is used for ≥18 m s−1 (Fairall et al. 2003). This might overestimate the roughness at very high wind speeds where Taylor and Yelland (2001) predict a decrease in the Charnock parameter.

With the availability of the improved COARE algorithm v3.0 (Fairall et al. 2003), it is now possible to obtain more realistic CD and CL values. This would reduce the uncertainties in the datasets (≈10%–20%) used for deriving the exchange coefficients presented in Kara et al. (2000, 2002). In addition, the use of the COARE v3.0 has the advantage of providing reliable transfer coefficients at low and high wind speeds. In this paper, we also provide an accurate formulation of air–sea fluxes using stability-dependent exchange coefficients, allowing the use of ocean model sea surface temperature (SST) in their calculation for a wide range of atmospheric conditions over the global ocean. This approach can substantially improve the accuracy of an ocean model simulated–SSTs by providing a physically natural tendency toward accurate SSTs without relaxation to any SST data or climatology.

The main purposes of this study are to present (i) exchange coefficients based on the COARE v3.0 (section 2) and (ii) polynomial parameterizations of CD, CL, and CS based on these coefficients (section 3). In particular, the stability correction in these parameterizations is approximated by a term depending on the (TaTs) and a polynomial function of the Va at 10 m using tabulated exchange coefficients as obtained from the COARE algorithm (v3.0). A virtual temperature-based correction to (TaTs) is later applied to the polynomial functions, thus including the effects of water vapor flux on the exchange coefficients (section 4). Finally, exchange coefficients presented in sections 3 and 4 are used for calculating climatological mean wind stress magnitude in the North Pacific Ocean over 1979–93 (section 5).

2. The COARE exchange coefficients

The COARE algorithm (v3.0) presented in Fairall et al. (2003) is based on previously published results and 2777 one-hour covariance flux measurements. Fairall et al. (2003) added 4439 new values from field experiments between 1997 and 1999 to test the algorithm. The new observational dataset now dominates the database, especially in the wind speed regime above 10 m s−1, where the number of observations increased from 67 to about 800. After applying quality controls, the database was used to evaluate the algorithm. The average (mean and median) model results agreed with the measurements to within about 5% for moisture from 0 to 20 m s−1. For stress, the covariance measurements were about 10% higher than the model at wind speeds over 15 m s−1 while inertial-dissipation measurements agreed closely at all wind speeds. Thus, it should be emphasized that the COARE algorithm (v3.0) is considered to be the state-of-the-art algorithm in this paper.

For the purpose of this paper, the COARE algorithm (v3.0) was modified to produce CD and CL at various (TaTs) and Va intervals over the global ocean. In particular, the algorithm was run using various values of Ts, relative humidity (RH) just above the sea surface, Ta, and Va. There is no cool skin or warm-layer calculation; thus, a skin sea surface temperature value is used in the algorithm. The use of a bulk value is an approximation, which allows one to perform the calculation without specifying radiative parameters or a time history of the fluxes. The roughness length for temperature is assumed equal to that for humidity in the algorithm; thus, CL is equal to CS in this paper.

In addition to (TaTs) and Va, the effect of any water vapor flux through RH must be taken into account in calculating the exchange coefficients at low winds, a topic which will be discussed in section 4 in detail. In this section the exchange coefficients are obtained for a relatively low Ts (10°C) and for humidity close to or at saturation (RH = 100%) to investigate the impact of only (TaTs) on CD and CL. Note that dividing the mixing ratio (qa) by the saturation mixing ratio (qsat) and then multiplying this number by 100 determines the RH (in %). A Ts value of 10°C was chosen for consistency with the COARE algorithm testing, and because almost no sensitivity to the choice of Ts was found. This is true since the only effect of the temperature in the COARE algorithm enters through the nonlinearity of the saturation vapor pressure for water vapor, which has some effect at 30°C and upward.

Exchange coefficients for CD and CL are first generated using an RH value of 100% for various (TaTs) ranging from −8° to 7°C: 7°, 6°, 5°, 4°, 3°, 2°, 1°, 0.75°, 0.5°, 0.25°, 0, −0.098°, −0.25°, −0.75°, −1°, −2°, −3°, −4°, −6°, −8°C. These values were chosen to represent a wide variety of (TaTs) intervals over the global ocean. It is noted that a (TaTs) value of ≈−0.098° gives the neutral stability value. The resolution of (TaTs) intervals was increased near the neutral stability value to better examine behavior of exchange coefficients with respect to stability. In addition, the 10-m wind speed values of 1–40 m s−1 with increments of 1 m s−1 were used to determine coefficients at both low and high wind speed conditions (e.g., tropical cyclones and hurricanes).

Figure 1 shows CD and CL values obtained from the COARE algorithm when RH is set to 100%. Both CD and CL experience large changes with wind speed at 10 m above sea surface, and they are also quite different with respect to (TaTs) values for a given wind speed. In general, CD is always greater than 2.0 × 10−3 at Va > 20 m s−1 for any given (TaTs) value (Fig. 1a). On the other hand, CL does not have any significant variability at high wind speeds (Fig. 1b). In particular, CL is almost constant with values ranging between ≈1.2 × 10−3 and ≈1.3 × 10−3 for any given (TaTs) values above 20 m s−1. There are strong variations in the exchange coefficients with respect to (TaTs) at very low wind speeds (Va < 3 m s−1), and this is especially evident at Va values of 1 and 2 m s−1 (Figs. 1c,d). The Va = 1 m s−1 cases are intended to show how variable CD and CL can be at very low wind speed regimes. Even the CD and CL values for Va = 1 m s−1 are quite different from those for Va = 2 m s−1. We actually obtained CD and CL values from the COARE algorithm at each 1 m s−1 interval up to 40 m s−1, but they are only shown at 2 m s−1 intervals in the figure for visual clarity. Although exchange coefficients can exist for Va = 0 m s−1, they are ignored in the analyses because air–sea fluxes are quite small under calm wind conditions.

Fig. 1.

Variations of exchange coefficients with respect to (TaTs) and Va values obtained from the COARE (v3.0) algorithm: (a) CD values plotted against the observed 10 m Va; (b) CL values plotted against the observed 10 m Va; (c) CD values plotted against the observed (TaTs), where Ta is at 10 m; and (d) CL values plotted against the observed (TaTs), where Ta is at 10 m. The small panels inside (a) and (b) are intended to show variability of CD and CL for only Va ≤ 10 m s−1, respectively. Dashed lines in (c) and (d) show CD and CL values for Va = 1 m s−1. It is noted that the bulk formulas are QL = CL L ρa Va (qaqs) for latent heat flux, Qs = CS Cp ρa Va (TaTs) for sensible heat flux, and τ = ρaCDV2a for wind stress. Here, L is the latent heat of vaporization (2.5 × 106 J kg−1), ρa = 100 Pa/[Rgas (Ta + 273.16)] (in kg m−3), where Rgas is the gas constant (287.1 J kg−1 K−1), and Pa is the mean sea level pressure.

Fig. 1.

Variations of exchange coefficients with respect to (TaTs) and Va values obtained from the COARE (v3.0) algorithm: (a) CD values plotted against the observed 10 m Va; (b) CL values plotted against the observed 10 m Va; (c) CD values plotted against the observed (TaTs), where Ta is at 10 m; and (d) CL values plotted against the observed (TaTs), where Ta is at 10 m. The small panels inside (a) and (b) are intended to show variability of CD and CL for only Va ≤ 10 m s−1, respectively. Dashed lines in (c) and (d) show CD and CL values for Va = 1 m s−1. It is noted that the bulk formulas are QL = CL L ρa Va (qaqs) for latent heat flux, Qs = CS Cp ρa Va (TaTs) for sensible heat flux, and τ = ρaCDV2a for wind stress. Here, L is the latent heat of vaporization (2.5 × 106 J kg−1), ρa = 100 Pa/[Rgas (Ta + 273.16)] (in kg m−3), where Rgas is the gas constant (287.1 J kg−1 K−1), and Pa is the mean sea level pressure.

The impact of stability through (TaTs) on both CD and CL is clearly evident from Fig. 1. Thus, exchange coefficients reduced to neutral stability, based solely on wind speed at 10 m above the sea surface (e.g., Large and Pond 1981), with no stability dependence (that is, TaTs = 0 and/or at saturated conditions), will not be appropriate when they are used in OGCM simulations. This becomes a critical issue in model simulations, especially when studying the diurnal cycle of upper ocean quantities.

Regarding hurricane-force winds (e.g., Va > 30 m s−1), there is almost no direct experimental evidence for determining exchange coefficients. There are very few direct observations of exchange coefficients for wind speeds greater than 25 m s−1. As to CL and CS, there is the problem of spray (e.g., Andreas and Emanuel 2001; Emanuel 2003). The COARE algorithm used here includes an option to define roughness length in terms of the sea state. It includes a formula based on wave age (Oost et al. 2002) or one based on significant height and slope of the waves (Taylor and Yelland 2001). The latter formula predicts significantly lower roughness at high wind speeds compared to the former one as explained by Taylor and Yelland (2003) who evaluated performance of the different formulas in predicting the CD. Based on evidence from laboratory studies (e.g., Keller et al. 1992) and comparisons with field experiments (e.g., Janssen 1997; Johnson et al. 1998), the formula proposed by Taylor and Yelland (2001), as used in the COARE algorithm, is applicable over most conditions. It also succeeds at high winds and very short fetches in wind-wave flumes when wave age formulas fail.

It should be emphasized that there is very little reason for the exchange coefficients obtained from the algorithm to change geographically. The main change from one area to another would be due to different stability or different temperature levels, but any such changes should be dealt with by the algorithm. The only other possible changes are 1) the change in gravity, which increases by up to 0.5% going from equator to the Pole; and 2) the nominal atmospheric boundary layer depth (600 m) used, which is also realistic in most areas. Another possible effect with geographic variation is the sea state. While the effect of sea state on surface roughness is debatable, it is usually negligible, say for Va > 10 m s−1 and in deep water, as mentioned in Taylor and Yelland (2001).

3. Polynomial approach for exchange coefficients

The exchange coefficients (CD and CL) are dependent on both Va and (TaTs) as discussed in section 2 (see Fig. 1). The effect of RH on CD and CL will be added to the polynomial functions as a correction term because including the full stability conditions (i.e., both TaTs and RH) in polynomial equations is not straightforward (see section 4).

We first investigate whether or not CD and CL can be expressed as simple polynomial functions of Va and (TaTs) using CD and CL values when humidity is close to or at saturation (i.e., RH = 100%). For this purpose, CD and CL values are plotted against (TaTs) for each Va value. As an example, Fig. 2 shows variation of exchange coefficients with respect to (TaTs) for Va = 1 m s−1 and Va = 6 m s−1, suggesting that we can construct polynomial functions in (TaTs) with coefficients that vary with Va. We split the whole (TaTs) range (from −8.00° to 7.00°C) into three subranges. The reason for this split was that there was no simple (low order) polynomial fit for CD and CL at low wind speeds when using the whole (TaTs) range.

Fig. 2.

Examples of applying optimal polynomial fits to the exchange coefficients: (a) CD for Va = 1 m s−1, (b) CL for Va = 1 m s−1, (c) CD for Va = 6 m s−1, and (d) CL for Va = 6 m s−1. Note that the exchange coefficients shown with open circles are obtained from the COARE algorithm.

Fig. 2.

Examples of applying optimal polynomial fits to the exchange coefficients: (a) CD for Va = 1 m s−1, (b) CL for Va = 1 m s−1, (c) CD for Va = 6 m s−1, and (d) CL for Va = 6 m s−1. Note that the exchange coefficients shown with open circles are obtained from the COARE algorithm.

Various (TaTs) intervals were tested to obtain the best available lowest-order polynomial function for the exchange coefficients. After a few tests, we found three ranges for (TaTs): 1) −8.00°C ≤ TaTs < −0.75°C, 2) −0.75°C ≤ TaTs ≤ 0.75°C, and 3) 0.75°C < TaTs ≤ 7.00°C. Hereinafter, for convenience, these intervals will be referred to as unstable, neutral, and stable cases, respectively. As evident from Fig. 2, it was possible to represent the neutral case using a linear equation, and the other two ranges (stable and unstable cases) by a quadratic equation for CD and CL. Obviously, a higher polynomial fit would fit better for the neutral case (e.g., Figs. 2a,b) but the use of the linear function is preferred for computational efficiency. The resulting error with respect to the COARE algorithm would also be negligible because a wind speed of 1 m s−1 entering the bulk formulations of wind stress and latent heat flux is also small for these two distinct cases.

We obtained a single equation for CD and CL that represents all wind speeds from Va values of 1–40 m s−1. Given that a quadratic polynomial is suitable for unstable and stable cases, and a linear fit is suitable for the neutral case, we then write a generalized polynomial equation for CD as follows:

 
formula

The same type of equations are also written for CL as follows:

 
formula

Unlike Kara et al. (2000, 2002), there are no limits set on Va in these formulations because we were able to take low and high wind speed values into account in the algorithms by using Va of 1, 2, and 3 m s−1 as obtained from the COARE algorithm.

The next step is to obtain CD0, CD1, and CD2 values (and similarly CL0, CL1, and CL2 values) to be used in the generalized equations of (1) and (2). To accomplish this, the example of the curve-fitting process shown in Fig. 2 for Va values of 1 and 6 m s−1 is repeated for each Va value ranging from 1 to 40 m s−1 (with 1 m s−1 increments). Thus, for example, the quadratic fits to CD for the unstable case yielded a total of 40 CD0, 40 CD1, and 40 CD2 values. Applying a linear fit would be sufficient for unstable and stable cases when Va ≥ 10 m s−1, but for consistency a quadratic fit was used for all Va values.

The final step is to express CD0, CD1, and CD2 values (similarly for CL0, CL1, and CL2 values) as polynomial functions of Va or V−1a, representing all Va values from 1 to 40 m s−1. However, there was no single low-order polynomial equation for the whole Va range using all 40 of the CD0, CD1, and CD2 values. This is not suprising because both CD and CL vary significantly at low wind speeds as shown before (see Figs. 1c,d). Therefore, we used two Va intervals: 1) Va ≤ 5 m s−1 and 2) Va ≥ 5 m s−1. Figure 3 shows CD0 and CL0 values for the Va ≤ 5 m s−1 case. While a cubic fit worked best for CD0, a quadratic fit is sufficient for CL0. In the unstable case, we can then write CD0 = 1.89100 − 0.71820 Va + 0.19750 V2a − 0.01790 V3a and CL0 = 2.07700 − 0.39330 Va + 0.03971 V2a.

Fig. 3.

The wind speed dependence values of CD0 and CL0. Results are shown for Va ≤ 5 m s−1 cases for each given (TaTs) interval.

Fig. 3.

The wind speed dependence values of CD0 and CL0. Results are shown for Va ≤ 5 m s−1 cases for each given (TaTs) interval.

Polynomial coefficients for CD and CL are given in Table 1 for unstable, neutral, and stable cases. Note that the best fit is sometimes obtained when using V−1a rather than Va as the dependent variable. Since V−1a is being used in some cases, there must be a minimum for Va. Here the minimum Va considered is 1 m s−1. It is also important to note that we obtained all polynomial functions such that they match up at Va = 5 m s−1, and at (TaTs) values of −0.75° and 0.75°C, as well.

Table 1.

Optimal polynomial coefficients for CD and CL. They are shown for Va ≤ 5 m s−1 and Va ≥ 5 m s−1 separately. Note that CD and CL coefficients are continuous at Va = 5 m s−1, (TaTs) = −0.75°C, and (TaTs) = 0.75°C.

Optimal polynomial coefficients for CD and CL. They are shown for Va ≤ 5 m s−1 and Va ≥ 5 m s−1 separately. Note that CD and CL coefficients are continuous at Va = 5 m s−1, (Ta − Ts) = −0.75°C, and (Ta − Ts) = 0.75°C.
Optimal polynomial coefficients for CD and CL. They are shown for Va ≤ 5 m s−1 and Va ≥ 5 m s−1 separately. Note that CD and CL coefficients are continuous at Va = 5 m s−1, (Ta − Ts) = −0.75°C, and (Ta − Ts) = 0.75°C.

4. Relative humidity effects on exchange coefficients

The polynomial functions presented in section 3 were determined using exchange coefficients from the COARE v3.0 algorithm, which was run for various (TaTs) values at RH = 100%, ignoring the impact of RH on the exchange coefficients. This simply implies that in our approximations so far, by far the major effect of water vapor on CD and CL is its temperature-dependent impact on stability, which is very large in the Tropics and decreases with decreasing absolute temperature. Because of the exponential increase of saturation vapor pressure with temperature, in regions such as the Tropics, humidity has a first-order effect on the stability, and since the Tropics tend to be light wind areas, on the fluxes as well.

Here, our purpose is to include full atmospheric stability in the polynomial functions (see section 3) by using water vapor in addition to (TaTs) and Va. The exchange coefficients calculated using these full stability conditions will be correct over a wide range of conditions. For example, our polynomial approach is applicable to conditions where the water vapor (including RH) makes a significant contribution to the buoyancy flux, that is, when the full impact of stability on the air–sea fluxes is included. Similarly, one would not expect any significant error when Ta is warmer than Ts, but the water vapor results in unstable stratification.

The polynomial functions for CD and CL shown in section 3 assume saturated conditions (i.e., RH = 100%). However, water vapor effects can be taken into account via a virtual temperature-based correction to (TaTs). The virtual air temperature (Tυa) is calculated from the mixing ratio at 10 m above the sea surface (qa) using Tυa = Ta [1 + 0.61 qa(Ta)] with Ta in Kelvins (K). One possible approach would be recast the polynomial functions for CD and CL to be functions of (TυaTυs) rather than (TaTs). Since we already have polynomial functions for saturated conditions, we instead apply a correction that is zero when RH = 100%:

 
formula

where qsat is the saturation mixing ratio at Ta from a simplified version of the original formulation for saturated vapor-pressure presented by Lowe (1977), which has a computational advantage over Buck (1981).

The polynomial functions for CD and CL applied to (TaTs)c are obviously identical to the original approach or saturated conditions and are in good agreement with the COARE algorithm (v3.0) for all relative humidities (not shown). Once the correction is applied to the polynomial functions, exchange coefficients are calculated. It is emphasized that in the polynomial functions one could use very fine increments for (TaTs) rather than those used in the original COARE algorithm (see section 2) to calculate CD and CL because all polynomials are continuous for any given (TaTs), Va, and RH.

Using the tabulated exchange coefficients obtained from the polynomial functions, we now investigate the impact of RH (0%–100%) in calculating CD and CL. It is clear that assuming an RH value of 100% results in significant errors for Va = 1 m s−1 (Fig. 4). The impact of RH on CD and CL is clearly seen when the difference between Ta and Ts is small (e.g., −2°C ≤ TaTs < 2°C). Such a result clearly demonstrates that atmospheric stability can be predicted incorrectly when ignoring RH dependence (i.e., assuming humidity close to or at saturation) in calculating CD and CL at low wind speeds. Overall, the direct effect of decreasing RH is to generally increase exchange coefficients as seen from the scatterplots of CD and CL for Va = 1 m s−1 (Fig. 5). For example, if RH were 20% and one assumed an RH value of 100%, then CD would be underestimated with ≈6 times less than its actual value (0.176 × 10−3 instead of 0.762 × 10−3) for Va = 1 m s−1 and TaTs = 1.0°C (Table 2). Similarly, the underestimate of CL would be ≈7 times less than its actual value (0.220 × 10−3 instead of 1.626 × 10−3) for Va = 1 m s−1 and TaTs = 0.7°C (Table 2). On the other hand, it must be noted that the effect of RH on exchange coefficients is relatively small for high wind speeds as can be seen in CD and CL values for Va = 6 m s−1 (Table 3). While the water vapor effect is most important at low wind speed, the reader is cautioned that low wind speed implies small wind stress and small latent and sensible heat fluxes. If the wind stress or latent and sensible heat flux values are small to start with, even an extremely large change in their strength might still mean small values after the RH effect is taken into account. For example, latent heat flux at Va = 1 m s−1 is ≈1.9 W m−2 (≈11.8 W m−2) when RH = 100% (RH = 20%). In both estimates we use TaTs = 2°C, and qaqs = 3 g kg−1.

Fig. 4.

Exchange coefficients plotted against various (TaTs) values at Va = 1 m s−1 as determined from the polynomial functions: (a) CD and (b) CL. This figure demonstrates variability of exchange coefficients with water vapor. Results are shown for RH values of 20%, 40%, 60%, 80%, and 100%.

Fig. 4.

Exchange coefficients plotted against various (TaTs) values at Va = 1 m s−1 as determined from the polynomial functions: (a) CD and (b) CL. This figure demonstrates variability of exchange coefficients with water vapor. Results are shown for RH values of 20%, 40%, 60%, 80%, and 100%.

Fig. 5.

Scatterplots of exchange coefficients for RH = 100% vs those for varying RH values at Va = 1 m s−1 (see also Fig. 4): (a) CD and (b) CL. The diagonal line identifies almost universal increase in exchange coefficients for given RH values against those for RH = 100%.

Fig. 5.

Scatterplots of exchange coefficients for RH = 100% vs those for varying RH values at Va = 1 m s−1 (see also Fig. 4): (a) CD and (b) CL. The diagonal line identifies almost universal increase in exchange coefficients for given RH values against those for RH = 100%.

Table 2.

Exchange coefficients for varying RH values and (TaTs) intervals at Va = 1 m s−1 as determined from the polynomial functions. Also given are CD (RH)/CD (RH = 100%) ratios. In the table, for example, CD (80%)/CD (RH = 100%) is the ratio of CD for RH = 80% to that for RH = 100%. Note that all exchange coefficients must be multiplied by 10−3.

Exchange coefficients for varying RH values and (Ta − Ts) intervals at Va = 1 m s−1 as determined from the polynomial functions. Also given are CD (RH)/CD (RH = 100%) ratios. In the table, for example, CD (80%)/CD (RH = 100%) is the ratio of CD for RH = 80% to that for RH = 100%. Note that all exchange coefficients must be multiplied by 10−3.
Exchange coefficients for varying RH values and (Ta − Ts) intervals at Va = 1 m s−1 as determined from the polynomial functions. Also given are CD (RH)/CD (RH = 100%) ratios. In the table, for example, CD (80%)/CD (RH = 100%) is the ratio of CD for RH = 80% to that for RH = 100%. Note that all exchange coefficients must be multiplied by 10−3.
Table 3.

Exchange coefficients for varying RH values and (TaTs) intervals at Va = 6 m s−1 as determined from the polynomial functions. Also given are CD (RH)/CD (RH = 100%) ratios. In the table, for example, CD (80%)/CD (RH = 100%) is the ratio of CD for RH = 80% to that for RH = 100%. Note that all exchange coefficients must be multiplied by 10−3.

Exchange coefficients for varying RH values and (Ta − Ts) intervals at Va = 6 m s−1 as determined from the polynomial functions. Also given are CD (RH)/CD (RH = 100%) ratios. In the table, for example, CD (80%)/CD (RH = 100%) is the ratio of CD for RH = 80% to that for RH = 100%. Note that all exchange coefficients must be multiplied by 10−3.
Exchange coefficients for varying RH values and (Ta − Ts) intervals at Va = 6 m s−1 as determined from the polynomial functions. Also given are CD (RH)/CD (RH = 100%) ratios. In the table, for example, CD (80%)/CD (RH = 100%) is the ratio of CD for RH = 80% to that for RH = 100%. Note that all exchange coefficients must be multiplied by 10−3.

5. Wind stress magnitude over the North Pacific Ocean

As an example, results from calculating wind stress magnitude based on the stability dependent CD is demonstrated in the North Pacific Ocean (Fig. 6). All analyses are based on 6-hourly meteorological data from an archived product, the 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40; Simmons and Gibson 2000). A climatological mean wind stress magnitude is formed over 1979–93. This time period was chosen to be consistent with ERA-15, which has been extensively used in previous ocean modeling studies (e.g., Hurlburt and Hogan 2000; Kara et al. 2004). In addition, all fields shown are interpolated onto the 1/12° Hybrid Coordinate Ocean Model (HYCOM) domain to demonstrate possible problems when using atmospheric variables (e.g., wind stress fields) from coarse-resolution archived products (e.g., 1.125°ERA-40 data) in forcing much finer OGCMs (e.g., 1/12° HYCOM). Basic features and a description of HYCOM are given in Kara et al. (2005).

Fig. 6.

Climatological mean wind stress magnitude over 1979–93 (a) as obtained directly from the ERA-40 output itself, (b) as computed from the bulk formulation using CD ignoring RH effects (see section 3), and (c) as computed from the bulk formulation using CD with the inclusion of RH effects (see section 4). Also shown are (d) wind speed at 10 m above the sea surface calculated from wind stress field shown in (c), (e) the ratio of wind stress magnitude calculated using (c)/(b), and (f) the land–sea mask from ERA-40, and note that a grid point is either land (1) or sea (0), giving values between 0 and 1 for interpolated the land–sea mask. All fields shown are interpolated on the 1/12° Mercator grid used in the Hybrid Coordinate Ocean Model (HYCOM).

Fig. 6.

Climatological mean wind stress magnitude over 1979–93 (a) as obtained directly from the ERA-40 output itself, (b) as computed from the bulk formulation using CD ignoring RH effects (see section 3), and (c) as computed from the bulk formulation using CD with the inclusion of RH effects (see section 4). Also shown are (d) wind speed at 10 m above the sea surface calculated from wind stress field shown in (c), (e) the ratio of wind stress magnitude calculated using (c)/(b), and (f) the land–sea mask from ERA-40, and note that a grid point is either land (1) or sea (0), giving values between 0 and 1 for interpolated the land–sea mask. All fields shown are interpolated on the 1/12° Mercator grid used in the Hybrid Coordinate Ocean Model (HYCOM).

We first form a climatological mean total wind stress magnitude field (Fig. 6a) using zonal and meridional wind stress components that are directly output at 6-hourly time intervals from ERA-40. This field is obtained from the planetary boundary layer model in the ECMWF atmospheric model. There is excessive noise in the wind stress magnitude near the coastal boundaries, such as along the eastern Pacific and in the Sea of Japan. This noise is a major motivation for recalculating the wind stress magnitude rather than using the one output from ERA-40. Thus, our purpose is to eliminate the noise shown in Fig. 6a, so that this field can be suitable as an atmospheric forcing field for an OGCM simulation. We recompute the wind stresses [τ = (τx, τy)] with a bulk formulation using CD values as derived in sections 3 and 4 as follows:

 
formula
 
formula
 
formula

where Va = (u2 + υ2)1/2, τ = (τ 2x + τ 2y)1/2, and u (υ) is the zonal (meridional) wind speed component (m s−1). Density of the air at the air–sea interface is obtained using ρa = 100 Pa/[Rgas (Ta + 273.16)] (in kg m−3), where Rgas is the gas constant (287.1 J kg−1 K−1) and Pa is the mean sea level pressure. All parameters needed for calculation of CD, ρa, and wind stress magnitude (i.e., Ta, Ts, qa, qsat, u, υ, Pa) are obtained from ERA-40 as mentioned earlier. Wind stress magnitude is computed at each 6-h interval and averaged over 1979–93.

The unrealistic noise along the coastlines is greatly cleaned up by using the bulk formulation with stability-dependent CD based solely on TaTs (Fig. 6b). This is also true when calculating CD depending on both TaTs and water vapor (Fig. 6c). The spatial patterns and approximate magnitude of the mean wind stress magnitude from the bulk formulation usually agree with those from ERA-40 itself. Basin-averaged mean wind stress magnitudes are 0.058, 0.034, and 0.037 N m−2 for Fig. 6a, Fig. 6b, and Fig. 6c, respectively. The relatively large basin-averaged mean value of 0.058 N m−2 in Fig. 6a is almost twice as large as those in Figs. 6b and 6c. This is due generally to the unrealistically large wind stress magnitudes near the land–sea boundaries, resulting from excessive noise in the original ERA-40 product.

The recomputed wind stress magnitude can also be used for determining Va. Thus, excessive noise near the coastal boundaries can be removed from the Va field as well (Fig. 6d). Here, assuming neutral CD values, Va is approximated using Va = 1/[ρa × 10−3 (0.926 + 3.218τ − 5.231τ 2 + 3.236τ3)] for τ ≤ 0.7711 and Va = 1/[ρa × 10−3 (1.461 + 0.485τ − 0.092τ 2 + 0.007τ3)] for τ > 0.7711. The use of the wind stress magnitude field in Fig. 6c instead of that in Fig. 6b did not make any significant difference in computing the Va because the ratio of these two fields is close to 1 over the most of North Pacific Ocean (Fig. 6e), revealing that the impact of water vapor in calculating CD is small on the climatological time scales. However, there are differences when Va is small, and this is consistent with the results shown in section 4. The reader should be cautioned that there could be impacts of water vapor on CD on shorter time scales (e.g., daily or shorter), and this is especially true at low wind speeds and low relative humidities (see section 4).

Finally, Fig. 6f shows the land–sea mask used in the original 1.125° × 1.125° Gaussian grid for ERA-40 interpolated onto the 1/12° HYCOM domain. For example, a contour value of 0.6 implies that the wind stress magnitudes on the model grid are ≈60% contaminated by the wind stress magnitude over land. Some grids near the coast are represented as land points in ERA-40 but most of them should be ocean. Such problems in atmospheric forcing fields were also noted other studies (e.g., Metzger and Hurlburt 2001; Kara et al. 2005). Thus, it is clear that unrealistic wind stress magnitudes near such coastal region errors (see Fig. 6a) are largely due to the land–sea mask and orographically generated Gibb’s waves. Gibb’s waves are caused by the spectral truncation of the orography to triangular truncation (P. Kållberg, ECMWF, 2003, personal communication). When averaging fields originally created in the ERA-40’s spectral geometry, traces of the originally minute Gibb’s waves become more visible. For example, Metzger (2003) noted such atmospheric forcing problems in an ocean model study of the South China Sea. On the other hand, the atmospheric variables (e.g., Ta, u, υ, etc) in the ERA-40 product are from postprocessed fields where the spectral traces have been removed, allowing one to calculate, for example, wind stress magnitude from these variables as we did in the North Pacific Ocean (Figs. 6b,c).

6. Summary and conclusions

Ocean models or coupled atmosphere–ocean models can be directly forced with net surface heat fluxes obtained from archived weather products. However, recalculating air–sea fluxes that are dependent upon the model SST is typically sufficient to keep model-simulated SSTs on track without any explicit relaxation to SST data or a climatology (although assimilation of SST data can further increase the model SST accuracy). This process requires reliable and efficient parameterizations of stability-dependent exchange coefficients over the global ocean, a topic that is investigated in this study.

We present simple and computationally efficient parameterizations of exchange coefficients for wind stress (CD) and latent heat flux (CL), which are of particular value not only for high-resolution global OGCM applications but also for various ocean mixed-layer and air–sea interaction studies. The polynomial functions for exchange coefficients are determined after running the COARE v3.0 algorithm for wind speed magnitudes at 10 m above the sea surface ranging from 1 to 40 m s−1. The stability dependence in the exchange coefficients was accomplished by using air–sea temperature difference and air–sea mixing ratio difference through RH for a given wind speed at 10 m above the sea surface. It is shown that ignoring water vapor effects in the calculation of exchange coefficients can result in serious errors at low wind speeds, generally a serious underestimation of the exchange coefficients for Va < 5 m s−1. As an application of our methodology, climatological mean wind stress magnitude is calculated using 6-hourly meteorological variables from ERA-40 in the North Pacific Ocean over 1979–93. Such an analysis can identify possible problems existing in atmospheric forcing parameters (e.g., wind stress magnitude), which are directly output from coarse-resolution archived weather products (e.g., ERA-40) and demonstrates how these problems can be reduced using our methodology. The use of exchange coefficients from the polynomial functions in the wind stress bulk formulation not only reproduces wind stress magnitudes approximately comparable to the original ERA-40 field in the interior of the North Pacific Ocean, but also eliminates unrealistic and excessive noise near the land–sea boundaries that exist in the original ERA-40 field.

Finally, the motivation for this paper was to present sufficiently accurate transfer coefficients for air–sea fluxes that were also computationally efficient enough to be used for climate model applications. Although the earlier COARE algorithm (v2.5b) required ≈20 iterations, the latest version (v3.0) iterates 3 times or less. However, even having three iterations at each time step would result in significant computational inefficiency for a fine resolution global OGCM that runs many years (e.g., decadal to interdecadal time scales). Exchange coefficients based on our derived polynomial functions can be calculated at ≈5 times lower computational cost than the more involved computations of COARE algorithm (v3.0) with three iterations for atmospheric stability. We also note here that the sensible heat flux coefficient CS is assumed to be equal to CL in this paper because, for simplicity, the roughness length for temperature is assumed equal to that for humidity in the algorithm. Further investigation concerning variable roughness lengths is beyond the scope of this study.

The exchange coefficients for air–sea fluxes presented in this paper are called the Naval Research Laboratory (NRL) Air–Sea Exchange Coefficients (NASEC). Two simple FORTRAN programs (one for CD and one for CL) along with the tabulated exchange coefficients based on various wind speed, relative humidity and air–sea temperature difference intervals can be downloaded from the Web site (available online at http://www7320.nrlssc.navy.mil/nasec/nasec.html).

Acknowledgments

The authors would like to extend special thanks to P. K. Taylor, Head of the James Rennell Division of the Southampton Oceanography Centre (SOC), for providing constructive criticisms to this paper. We appreciate the helpful comments of E. J. Metzger (NRL) during the preparation of this manuscript. We acknowledge the excellent comments of the two anonymous reviewers. This work is a contribution to the 6.2 project Hybrid Coordinate Ocean Model (HYCOM) and advanced data assimilation sponsored by the Office of Naval Research (ONR) under Program Element 602435N and to the project HYCOM consortium for data-assimilative ocean modeling sponsored under the National Ocean Partnership Program.

REFERENCES

REFERENCES
Andreas
,
E. L.
, and
K A.
Emanuel
,
2001
:
Effects of sea spray on tropical cyclone intensity.
J. Atmos. Sci.
,
58
,
3741
3751
.
Buck
,
A L.
,
1981
:
New equations for computing vapor pressure and enhancement factor.
J. Appl. Meteor.
,
20
,
1527
1532
.
Bunker
,
A F.
,
1976
:
Computations of surface energy flux and annual air–sea interaction cycles of the North Atlantic Ocean.
Mon. Wea. Rev.
,
104
,
1122
1140
.
DeCosmo
,
J.
,
K B.
Katsaros
,
S D.
Smith
,
R J.
Anderson
,
W A.
Oast
,
K.
Bumke
, and
H.
Chadwick
,
1996
:
Air–sea exchange of water vapor and sensible heat: The Humidity Exchange Over the Sea (HEXOS) results.
J. Geophys. Res.
,
101
,
12001
12016
.
Emanuel
,
K.
,
2003
:
A similarity hypothesis for air–sea exchange at extreme wind speeds.
J. Atmos. Sci.
,
60
,
1420
1428
.
Fairall
,
C W.
,
E F.
Bradley
,
D P.
Rogers
,
J B.
Edson
, and
G S.
Young
,
1996
:
Bulk parameterization of air–sea fluxes for Tropical Ocean–Global Atmosphere Coupled–Ocean Atmosphere Response Experiment.
J. Geophys. Res.
,
101
,
3747
3764
.
Fairall
,
C W.
,
E F.
Bradley
,
J E.
Hare
,
A A.
Grachev
, and
J B.
Edson
,
2003
:
Bulk parameterization of air–sea fluxes: Updates and verification for the COARE algorithm.
J. Climate
,
16
,
571
591
.
Hare
,
J E.
,
P. O. G.
Persson
,
C W.
Fairall
, and
J B.
Edson
,
1999
:
Behavior of Charnock’s relationship for high wind conditions. Preprints, 13th Symp. on Boundary Layers and Turbulence, Dallas, TX, Amer. Meteor. Soc., 252–255
.
Hurlburt
,
H E.
, and
P J.
Hogan
,
2000
:
Impact of 1/8° to 1/64° resolution on Gulf Stream model-data comparisons in basin-scale subtropical Atlantic Ocean models.
Dyn. Atmos. Oceans
,
32
,
283
329
.
Isemer
,
H-J.
,
J.
Willebrand
, and
L.
Hasse
,
1989
:
Fine adjustment of large scale air– sea energy flux parameterizations by direct estimates of ocean heat transport.
J. Climate
,
2
,
1173
1184
.
Janssen
,
J. A. M.
,
1997
:
Does wind stress depend on sea-state or not?—A statistical error analysis of HEXMAX data.
Bound.-Layer Meteor.
,
83
,
479
503
.
Johnson
,
H K.
,
J.
Hoejstrup
,
H J.
Vested
, and
S E.
Larsen
,
1998
:
Dependence of sea surface roughness on wind waves.
J. Phys. Oceanogr.
,
28
,
1702
1716
.
Kara
,
A B.
,
P A.
Rochford
, and
H E.
Hurlburt
,
2000
:
Efficient and accurate bulk parameterizations of air–sea fluxes for use in general circulation models.
J. Atmos. Oceanic Technol.
,
17
,
1421
1438
.
Kara
,
A B.
,
P A.
Rochford
, and
H E.
Hurlburt
,
2002
:
Air–sea flux estimates and the 1997–1998 ENSO event.
Bound.-Layer Meteor.
,
103
,
439
458
.
Kara
,
A B.
,
H E.
Hurlburt
,
P A.
Rochford
, and
J J.
O’Brien
,
2004
:
The impact of water turbidity on the interannual sea surface temperature simulations in a layered global ocean model.
J. Phys. Oceanogr.
,
34
,
345
359
.
Kara
,
A B.
,
A J.
Wallcraft
, and
H E.
Hurlburt
,
2005
:
Sea surface temperature sensitivity to water turbidity from simulations of the turbid Black Sea using HYCOM.
J. Phys. Oceanogr.
,
35
,
33
54
.
Keller
,
M R.
,
W C.
Keller
, and
W J.
Plant
,
1992
:
A wave tank study of the dependence of X band cross sections on wind speed and water temperature.
J. Geophys. Res.
,
97
,
5771
5792
.
Large
,
W G.
, and
S.
Pond
,
1981
:
Open ocean momentum flux measurements in moderate to strong winds.
J. Phys. Oceanogr.
,
11
,
324
336
.
Lowe
,
P R.
,
1977
:
An approximating polynomial for the computation of saturation vapor pressure.
J. Appl. Meteor.
,
16
,
100
103
.
Metzger
,
E J.
,
2003
:
Upper ocean sensitivity to wind forcing in the South China Sea.
J. Oceanogr.
,
59
,
783
798
.
Metzger
,
E J.
, and
H E.
Hurlburt
,
2001
:
The importance of high resolution and accurate coastline geometry in modeling South China Sea inflow.
Geophys. Res. Lett.
,
28
,
1059
1062
.
Oost
,
W A.
,
G J.
Komen
,
C. M. J.
Jacobs
, and
C.
van Oort
,
2002
:
New evidence for a relationship between wind stress and wave age from measurements during ASGAMAGE.
Bound.-Layer Meteor.
,
103
,
409
438
.
Simmons
,
A J.
, and
J K.
Gibson
,
2000
:
The ERA-40 project plan. ERA-40 Project Report Series, No. 1, 62 pp. [Available from ECMWF, Shinfield Park, Reading RG2 9AX, United Kingdom.]
.
Smith
,
S D.
,
1988
:
Coefficients for sea surface wind stress, heat flux and wind profiles as a function of wind speed and temperature.
J. Geophys. Res.
,
93
,
15467
15472
.
Taylor
,
P K.
, and
M J.
Yelland
,
2001
:
The dependence of sea surface roughness on the height and steepness of the waves.
J. Phys. Oceanogr.
,
31
,
572
590
.
Taylor
,
P K.
, and
M.
Yelland
,
2003
:
Sea surface roughness: High winds and short fetches. Preprints, 83d AMS Annual Meeting, Long Beach, CA, Amer. Meteor. Soc., CD-ROM, P3.5
.
Wallcraft
,
A J.
,
A B.
Kara
,
H E.
Hurlburt
, and
P A.
Rochford
,
2003
:
The NRL Layered Global Ocean Model (NLOM) with an embedded mixed layer submodel: Formulation and tuning.
J. Atmos. Oceanic Technol.
,
20
,
1601
1615
.
Yelland
,
M J.
,
B I.
Moat
,
P K.
Taylor
,
R W.
Pascal
,
J.
Hutchings
, and
V C.
Cornell
,
1998
:
Wind stress measurements from the open ocean corrected for airflow distortion by the ship.
J. Phys. Oceanogr.
,
28
,
1511
1526
.
Yelland
,
M J.
,
B I.
Moat
,
R W.
Pascal
, and
D I.
Berry
,
2002
:
CFD model estimates of the airflow distortion over research ships and the impact on momentum flux measurements.
J. Atmos. Oceanic Technol.
,
19
,
1477
1499
.

Footnotes

Corresponding author address: Birol Kara, Naval Research Laboratory, Code 7320, Bldg. 1009, Stennis Space Center, MS 39529-5004. Email: kara@nrlssc.navy.mil

* Naval Research Laboratory Contribution Number NRL/JA/7304/03/53.