JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 101, NO. C2, PAGES 3747-3764, FEBRUARY 15, 1996
Bulk parameterization of air-sea fluxes for Tropical Ocean-
Global Atmosphere Coupled-Ocean Atmosphere Response
Experiment
C. W. Fairall, • E. F. Bradley, • D. P. Rogers, 3 J. B. Edson, 4 and G. S. Young s
Abstract. This paper describes the various physical processes relating near-surface
atmospheric and oceanographic bulk variables; their relationship to the surface fluxes of
momentum, sensible heat, and latent heat; and their expression in a bulk flux algorithm.
The algorithm follows the standard Monin-Obukhov similarity approach for near-surface
meteorological measurements but includes separate models for the ocean's cool skin
and the diurnal warm layer, which are used to derive true skin temperature from the bulk
temperature measured at some depth near the surface. The basic structure is an out-
growth of the Liu-Katsaros-Businger [Liu et al., 1979] method, with modifications to
include a different specification of the roughness/stress relationship, a gustiness velocity
to account for the additional flux induced by boundary layer scale variability, and profile
functions obeying the convective limit. Additionally, we have considered the contribu-
tions of the sensible heat carried by precipitation and the requirement that the net dry
mass flux be zero (the so-called Webb correction [Webb et al., 1980]). The algorithm
has been tuned to fit measurements made on the R/V Moana Wave in the three different
cruise legs made during the Coupled Ocean-Atmosphere Response Experiment. These
measurements yielded 1622 fifty-min averages of fluxes and bulk variables in the wind
speed range from 0.5 to 10 rn s -•. The analysis gives statistically reliable values for the
Charnock [1955] constant (Ix = 0.011) and the gustiness parameter ([• = 1.25). An over-
all mean value for the latent heat flux, neutral bulk-transfer coefficient was 1.11 x 10 -3,
declining slightly with increasing wind speed. Mean values for the sensible and latent
heat fluxes were 9.1 and 103.5 W m-2; mean values for the Webb and rain heat fluxes
were 2.5 and 4.5 W m -2. Accounting for all factors, the net surface heat transfer to the
ocean was 17.9 + 10 W m -2.
1. Introduction
The importance of air-sea interaction to the Earth's
climate is widely appreciated. The pivotal role of the
tropical oceans in climate and interannual climate variability
led to the establishment of the Tropical Ocean-Global
Atmosphere (TOGA) program; the subsequent identification
of the dominance of the Pacific Ocean in this variability
resulted in the Coupled Ocean-Atmosphere Response
Experiment (COARE) [World Climate Research Program
(WCRP), 1990; Webster and Lukas, 1992]. The interfacial
fluxes are one of three elements of the COARE program.
The COARE science plan [WCRP, 1990, pp. A7-A8]
•Environmental Technology Laboratory, NOAA,
Boulder, Colorado.
2Centre for Environmental Mechanics, Commonwealth Scientific
and Industrial Research Organization, Canberra, Australian Capital
Territory, Australia.
3physical Oceanography Research Division, Scripps Institute
of Oceanography, La Jolla, California.
nApplied Ocean Physics and Engineering Department, Woods
Hole Oceanographic Institution, Woods Hole, Massachusetts.
5Department of Meteorology, Pennsylvania State University,
University Park. •
Copyright 1996 by the American Geophysical Union.
Paper number 95JC03205.
0148-0227/96/95JC-03205 $05.00
identifies several fundamental gaps in our knowledge
relevant to fluxes:
Atmospheric response in models is extremely
sensitive to SST variations, especially where SST is
warm. However, ocean models almost universally
predict temperatures which are too warm, probably
associated with poor assessments of heat, momentum,
and moisture fluxes to the ocean and atmosphere.
The heat balance of the warm pool region of the
western Pacific is poorly known with discrepancies as
large as 80 W m '2. The relative involvement of the
slowly evolving atmosphere or the higher-frequency,
more episodic, equatorial events is not understood.
Webster and Lukas [1992, p. 1394] emphasized that "the
variation of fluxes between the ocean and the atmosphere is
very sensitive to the choice of parameterization, especially in
low wind regimes." This has been verified by Miller et al.
[1992], who found dramatic improvements in simulated
tropical phenomena by strengthening the air-sea coupling in
the light-wind regime. Furthermore, observational problems
in the climatological database, particularly the air-sea
temperature difference [Lukas, 1989], the strong boundary
layer diurnal cycle in light winds, and the unknown sensible
cooling associated with precipitation, represent additional
uncertainties in assessing the surface energy balance of the
warm pool. In summary, our ability to diagnose, simulate,
3747
3748 FAIRALL ET AL.' BULK PARAMETERIZATION OF AIR-SEA FLUXES
and predict climate and climate variability is impaired by a
general lack of high-quality data in the region and inadequate
parameterizations of air-sea fluxes.
1.1. Scope of This Paper
This paper is concerned with the estimation of air-sea
fluxes from bulk variables, with a focus on a specific
algorithm developed for the TOGA COARE investigators.
Following background material on fluxes and similarity
theory, we discuss the representation of the near-surface
transfer processes in terms of the surface roughness para-
meters; theoretical issues associated with extending tradi-
tional methods to the light-wind, convective regime; the
proper thermodynamic constants for the computation of the
fluxes; and the estimation of the rainfall contribution to
surface cooling (section 2). In section 3 we touch briefly on
several measurement issues: the special problem of flux
measurements from ships, the accuracy requirements for bulk
variables, and a method to correct bulk water temperatures
to obtain the true interfacial sea surface temperature (SST).
The actual procedure used in the present algorithm is found
in section 4.2; readers familiar with the theory of turbulent
flux measurement might choose to refer immediately to this
section and use it as a guide to the preceding development.
A comparison of the algorithm with the Moana Wave
COARE data is given in section 5, including an analysis of
the relative contributions to the latent heat flux by stability
and gustiness corrections. Our conclusions are given in
section 6.
1.2. Bulk Flux Estimations
Simultaneous flux and bulk meteorological variable
measurements combined with laboratory studies of air-sea
transfer processes are used to develop the bulk formulas and
transfer coefficients. The classic reviews on this subject
[Garratt, 1977; Smith, 1988] reveal a substantial midlatitude
bias in the field measurements, with most of the data
obtained in the 4-15 m s '1 wind speed regime. The majority
of the data is also from offshore towers, coastal areas, and
other shallow water regimes. Uncertainties in the average
neutral coefficients as a function of wind speed remain the
critical question. Blanc's [1985] study showed a factor of 2
variation in suggested values for the humidity transfer
coefficient with a consensus uncertainty of about 30%. Note
that a 10% uncertainty in this transfer coefficient results in
a 10 W m '2 uncertainty in the latent heat flux and thus the
surface energy budget of the warm pool.
In low wind speed regimes it is necessary to account for
buoyancy effects on turbulent transport, and standard
stability-dependent bulk schemes [e.g., Liu et al., 1979;
Smith, 1988] have shown good performance in the tropics
[Bradley et al., 1991]. However, a careful analysis [Godfrey
and Beljaars, 1991] has shown these schemes to become
singular at winds speeds below about 0.5 m s 4. This occurs
when a basic similarity profile assumption (that the rough-
ness length is much smaller than the Monin-Obukhov length)
is violated. Godfrey and Beljaars [1991] showed that this
singularity can be eliminated by adding a "gustiness" velocity
wg related to the normal convective scaling velocity, which
accounts for the fact that the amplitude of the mean wind
vector does not properly characterize the mean wind speed
in light winds. Because 1-hour average point winds of less
than 4 m S -1 OCCUr about half of the time in the COARE
region, particular attention must be paid to this problem.
Recent experimental studies in the COARE region have
shown unequivocally that the scalar fluxes do not go to zero
in the limit of eero mean wind [Bradley et al., 1991;
Fujitani, 1992; Young et al., 1992; Bradley etal., 1993;
Greenhut and Khalsa, 1995].
2. Theory
2.1. Background
The turbulent fluxes of sensible heat H.•, latent heat H l,
and stress 'c components are defined by the normal Reynolds
averages,
H,,. = PaCpa w'T' = -PaCpa u, T,
H l = PaLe w'q' = -PaLeu, q.
(1)
2
• = Pa W'U' =-Pa u*
where w', T', q', and u ' represent the turbulent fluctuations
of vertical wind, temperature, water vapor mixing ratio, and
the streamwise component of horizontal wind, respectively;
T,, q,, and u, are the related Monin-Obukhov similarity
(MOS) scaling parameters [Panofsky and Dutton, 1984;
Geernaert, 1990]. The overbar denotes an ensemble average
but, in practice, is usually a time or space average.
The standard bulk expressions for the scalar fluxes and
stress components are
H.,. = PaCpa Ch S (T.• - O)
Hi = P a Le Ce S (q,,. - q)
'lj i -- PaCa S (usi- ui)
(2a)
(2b)
(2c)
where Ca, Ch, and C e are the transfer coefficients for stress,
sensible heat, and latent heat, respectively; 0 is the potential
temperature, q is the water vapor mixing ratio, and u• is one
of the horizontal wind components relative to the fixed
Earth, each measured at some atmospheric reference height
z r and averaged as in (1). S is the average value of the wind
speed relative to the sea surface at Zr; T.• is the sea surface
interface temperature; u.• i is the surface current; and q,• is the
interfacial value of the water vapor mixing ratio that is
computed from the saturation mixing ratio for pure water at
the SST,
q,. = 0.98 qsat (T•.) (3)
Alternatively, we may measure the wind components relative
to the sea surface, in which case the u.,. i terms are zero.
Following Sverdrup et al. [ 1942], the factor of 0.98 multiply-
ing the saturation specific humidity of the SST takes into
account the reduction in vapor pressure caused by a typical
salinity of 34 parts per thousand. Note that
0 = T + 0.0098 Z r
q = RHq.•(T)
(4)
where T is the air temperature at z r and RH is the relative
humidity.
FAIRALL ET AL.' BULK PARAMETERIZATION OF AIR-SEA FLUXES 3749
10
x
x
1 2 3 4 5 6 7 8 9
U (ms '1)
I I ,
lO 11 12
Figure 1. Roughness Reynolds number R r as a function of
the 10-m wind speed for the R/V Moana Wave Coupled
Ocean-Atmosphere Response Experiment (COARE) data.
The data have been averaged in wind speed bins with 1 m s 4
bin widths. The solid line is the COARE 2.0 algorithm
result, and the crosses are computed from the mean inertial-
dissipation values for stress.
The transfer coefficients
individual profile components,
in (2) are partitioned into
1/2 1/2
C h - CT C d
•/2 •/2 (5)
C e - Cq C d
1/2 1/2
C d ---- C d C d
which are themselves functions of the fluxes in a manner
described by MOS surface-layer theory [Panofsky and
Dutton, 1984; Geernaert, 1990],
/ Crn
C T = CTn
aK
1/2 1/2 / Cqn
C q = C qn I - •/h (•) (6)
aK
1/2 1/2 1 - •u(•)
C d = Cdn
Here < is the von Kfirmfin constant (0.4), a accounts for the
difference in scalar and velocity von Kfirmfin constants, • is
the MOS profile function (assumed the same for temperature
and humidity), and • = z r/L, where
L_ • = Kg
r (r. + 0.6 rq.)/u (7)
The subscript n denotes the value in neutral conditions
(i.e., • = 0) where • = 0. For a reference height of 10 m,
Sen and Chn are approximately 1 x 10 -3 and have little wind
speed dependence [Liu et al., 1979; Smith, 1988]. The
neutral transfer coefficients are related to the roughness
lengths (Zo for velocity, Zor for temperature, and Zoq for
humidity), which are defined as the height where the
extrapolation of the log-z portion of the respective profile
(of u, T, or q) intersects the surface value:
•/2 a K
CTn
10g(Zr/Zor)
•/2 a K
C qn 10g(Zr/Zoq) (8)
Cdn =
10g(Zr/Z o)
The scaling parameters from (1) can be computed
independently from the transfer coefficients given in (6):
1/2
T, = - c (r.,.- 0)
1/2 _
q, = -cq (q,. q)
2
u, = C a Su
(9)
where u denotes the magnitude of the mean wind vector
(relative to the sea surface).
2.2. Surface Characterization
The velocity roughness length Zo is often crudely related
to the physical roughness of the surface [see Panofsky and
Dutton, 1984, p. 123], but the scalar roughness lengths are
more complicated. This is discussed in detail by Garratt
[1992, chapters 4 and 5] or Kraus and Businger [1994,
chapter 5], so only a brief background will be given here.
From laboratory studies it has proven convenient to charac-
terize the surface and the flow regime by the roughness
Reynolds number,
U, Zo (10)
g r =
v
where v is the kinematic viscosity of air. For later use we
also define the scalar equivalents of gr: R r = (U, Zor)/V for
temperature and R• = (u, Zoo)Iv for humidity. Figure 1
shows the relationship between R r and wind speed, obtained
from Moana Wave data during COARE, which will be
described in detail in section 5.2. According to these
classical studies [e.g., Kraus and Businger, 1994], when
g r < 0.13, the flow is said to be "aerodynamically smooth";
that is, the actual roughness elements on the surface are
irrelevant and the surface stress is supported by viscous
shear. As the wind speed decreases, g r approaches a
constant value of about 0.11 and the relationship between
roughness and stress is fixed:
0.11v
z ø = (11)
For R r > 2.0 the flow is "rough" and the stress is dominated
by pressure and viscous transfers associated with the rough-
ness elements.
3750
FAIRALL ET AL.: BULK PARAMETERIZATION OF AIR-SEA FLUXES
Over the ocean, smooth flow occurs for 10-m wind speeds
less than about 2 m s '• and rough flow occurs for wind speed
greater than about 8 m s '•. The roughness elements over the
ocean are primarily surface gravity waves that are generated
by the win•stress. On the basis of scaling arguments about
the slope of the average locally generated seas as a function
of the surface stress, Charnock [1955] gave the mean
relationship between oceanic roughness and stress for rough
flow:
2
{•U,
Zo = (12)
g
where • is the "Charnock" constant for which values
between 0.010 and 0.035 can be found in the literature [e.g.,
Garratt, 1992, Table 4.1]. The value of the Charnock
"constant" can be linked to gross characterizations of the sea
state [Geernaert, 1990; Nordeng, 1991] such as the age or
slope of the dominant wavelength (from the peak of the
gravity wave spectrum). This has been used to explain, for
example, the increase of the velocity transfer coefficient in
shallow water. The literature on this subject is quite con-
flicting, so at this stage the subject must still be considered
exploratory. Broader application also awaits the ready
availability of wave spectral information. There have also
been suggestions [Wu, 1968] that capillary waves contribute
significantly to the stress at intermediate wind speeds. This
implies that the surface tension of the surface of the ocean
must be considered in wind/stress relationships. So far, clear
experimental verification is sketchy.
For 10-m wind speeds between 2 and 20 m s -j, decades of
field programs have not succeeded in clearly demonstrating
an open-ocean, neutral transfer coefficient for heat and
moisture that is significantly different than 1.1 x 10 '3 + 15%
[Garratt, 1992]. Despite this, the trend in characterizing heat
and moisture transfer has been to follow laboratory and
overland studies in parameterizing the scalar roughness
lengths in terms of the roughness Reynolds number. The
reasons are that the 15% uncertainty is no longer acceptable;
also, considerations of the roughness structure give us an
approach that can be extended to nonequilibrium wave states
and can also be used to deal with air-sea transfer of trace
gases. The leading examples are the model of Liu et al.
[ 1979] (hereinafter referred to as LKB) and Brutsaert [ 1982].
In the laboratory, simple experiments have been done to
determine the transfer of heat and moisture by molecular
diffusive processes in the thin sublayer directly adjacent to
the water surface. Brutsaert [1982] assumed that this near-
surface profile must match the log profile of turbulent
transfer at some matching height. This leads to a relation-
ship between the scalar roughness length (i.e., the log-profile
variable) and the velocity roughness length that depends on
Rr (equation (10)). The LKB model is quite similar to that
of Brutsaert, except the sublayer profile is given a specific
(exponential) shape and the matching height is determined
by equating the slopes of the two profile forms where they
intersect. This difference, plus alternative choices of
sublayer constants and specifications of the velocity
behavior, leads to substantial differences between the two
models with respect to the exchange coefficients used.
Figure 2 illustrates this with the Moana Wave COARE
measurements, which again will be described in detail in
section 5.2.
1.45
1.40
1.35
1.30
1.25
1.20
- \ /
\ /
ß
1.10 "..•••',." x x x
1.05 r I '1, I _1 I I I I I I I
1'00•) 1 2 3 4 5 6 7 $ 9 10 11 12
U (ms '1 )
Figure 2. Neutral stability values for the 10-m moisture
transfer coefficient Cen as a function of 10-m wind speed,
bin averaged as in Figure 1. The solid line is the
COARE 2.0 model, the dotted line is the Garratt/Brutsaert
model [Garratt, 1992; Brutsaert, 1982], and the dashed
line is the original Liu-Katsaros-Businger (LKB) model
[Liu et al., 1979] (see text). The crosses are derived from
the Moana Wave covariance latent heat flux measurements.
2.3. Convective Behavior
Within the framework of MOS many dynamical variables
have clearly defined asymptotic behavior in the so-called free
convection limit [Panofsky and Dutton, 1984; Garratt, 1992],
when u, approaches zero but the buoyancy flux does not
(that is, L goes to zero). For example, the dimensionless
vertical gradients of scalar quantities are expected to exhibit
a •-]/3 dependence. This leads to a scalar profile function of
the form
•. = 1.5In ,Y2+Y+I - •arctan +
3
(13)
where
y = •1- ¾• (14)
and ¾ is an empirical constant. Note that the convective limit
argument says nothing about the form of these functions near
neutral stability. Numerous overland field programs have
determined the forms of the profile functions for near-neutral
conditions [Hogstrom, 1988], but the experimental difficulties
of measuring the small gradients in the convective limit have
prevented clear verification of (13).
It is clear, however, that expressions such as (6) may
become singular if W becomes too large. Furthe.rmore,
Godfrey and Beljaars [ 1991 ] have pointed out that (6) cannot
be applied in the strict convective limit because it is based
on the constraint that the magnitude of Zoq/L << 1. This has
FAIRALL ET AL.: BULK PARAMETERIZATION OF AIR-SEA FLUXES 375!
resulted in attempts to sidestep the MOS framework and
scale the problem directly (see Liu [1989] for a discussion).
Laboratory studies have indicated that sensible heat flux does
scale as the 4/3 power of the air-sea temperature difference.
These studies use scaling arguments that depend on the
molecular diffusivities. These concepts have been general-
ized and applied in a limited way to the ocean [Golitsyn and
Grachev, 1986]. Following an approach based on boundary
layer convective similarity, Stull [1994] has developed a
scaling model where the sensible heat flux scales as the
3/2 power of temperature difference and the 1/2 power of the
depth of the convective boundary layer zt. Note that Stull's
theory was originally for smooth flow, where the actual
physical roughness of the surface is no longer relevant
(that is, zt is the length scale), but he has suggested that it is
also valid for rough flow.
Another approach has been developed by noting that in (2)
the parameter S is, in fact, the average value of the wind
speed, not the magnitude of the mean wind vector.
Schumann [1988] and Godfrey and Beljaars [1991]
expressed S as
2 2 2 u 2 2 (15)
= = + Wg
S 2 tt x + try + Wg
where u x and uy are the mean wind components, Wg is
proportional to the convective scaling velocity
Wg = [SW, (16)
and [5 is an empirical constant, of the order of 1.0, but which
depends on the temporal/spatial scale used to compute the
averages. We compute W, as follows:
W• g l H' HI I (17)
= " + 0.61 T z i
'• p,, c•,, p,,L•
Algebra shows that as the mean vector wind approaches
zero, 3' will approach [• W. and (2) will yield a result equiva-
lent to Stull's [ 1994] scaling theory.
Sykes et al. [1993] carried this concept further by examin-
ing the structure of the local profiles within the gusts.
Whereas (2c) properly implies that the average stress vector
approaches zero as the average vector wind approaches zero,
the local turbulence intensity that drives the scalar fluxes and
determines the local roughness length within the gusts must
be scaled by the wind speed:
to seawater), etc. Businger [1982] discussed at length the
issue of the additional contributions to the heat carried by
atmospheric moisture. He showed that the total enthalpy
transported by turbulent correlations is
p•,wh : p•, w'T' [cp•,+ q (c•v- c•,)]
+ p,w'q' [Ct•v(r-rr)+me]
(19)
where p, is the density of moist air, %, is the specific heat
of dry air, %v is the specific heat of water vapor, and Tr is
a reference temperature. This development follows Frank
and Emmitt [1981] and corrects an erroneous expression
developed by Brook [1978]. Businger [1982] pointed out
that the proper reference temperature in this application is the
SST, T,.. The second %• term represents the heat required to
cool the water vapor from T.,. to the air temperature after it is
evaporated. Thus the total sensible heat is
H.,. = pc•, {w'T' [1+ q(c•- c•,)/c•,]
+ w'q' (T- T.,) c•v/c• ' } (20)
For TOGA COARE both of these correction terms are, on
average, less than 0.2 W m '2 and of opposite sign. Because
it is unlikely that bulk fluxes will approach this accuracy,
even in an average sense, in the near future we suggest that
these correction terms be neglected.
A second issue is the proper lower boundary condition for
the total mass flux, as discussed by Webb et al. [1980],
which has become known as the "Webb effect." Webb et al.
[1980, p. 87] described the problem,
If the heat flux is upwards (positive) then rising air
parcels are on average warmer than descending parcels,
so that on the assumption of zero mean vertical mass
flow of air there must exist a small mean upward velo-
city component. Thus, in measurements which include
the fluctuations of density (of a minor constituent) and
of vertical airspeed, w, about its mean w'-, the contri-
bution to the flux of (the constituent) associated with
• is missed, and an appropriate correction having the
same sign as the heat flux must be added.
This small mean vertical velocity is given by
• = 1.61 w'q' + (1 + 1.61q) w'T'/r (21)
2 = Cd S 2 (18)
•,t
The point is if the wind blows at 1 m s -• from the east for
half an hour and then blows at 1 m s -1 from the west for half
an hour, it is the average wind speed of 1 m s -1 that must be
used to compute roughness Reynolds number, the scalar
transfer coefficients, and the MOS stability parameters.
This leads to a correction that must be added to the latent
heat flux:
Hlw = Pa Le •q (22)
For the COARE observation period, Hlw has an average
value of 4 W m -2.
2.4. Flux/Moisture Corrections
The flux accuracy guidelines set by the COARE working
groups have placed unprecedented demands on measurements
and computation. We no longer have the luxury of being
able to use "ballpark" figures for geophysical parameters
such as the acceleration of gravity, the latent heat of vapori-
zation of water, the vapor pressure of pure water (as opposed
2.5. Precipitation Effects
Gosnell et al. [1995] showed that the sensible heat
transferred to the ocean surface by the rain can be repre-
sented as
Hsr -' -gcpwOtw(l+ B; 1) AT (23)
where R is the rain rate (liquid water flux), %w is the specific