profile γx(d) is also required. This term arises from an analysis of the turbulent temperature-variance budget (Deardorff, 1972; Holtslag and Moeng, 1991) and from observations in the atmosphere that show non-zero buoyancy fluxes where the local buoyancy gradient is zero or even weakly of the opposite sign (i.e., a local deviation from the Fickian form (1) is required to avoid large negative µ values). Throughout the boundary layer, h ≥ d ≥ 0, the fluxes are given by
The diffusivity profile is formulated to conform to known and desirable properties of boundary-layer turbulence. In the surface layer, (with ), the semi-empirical results of Monin-Obukhov similarity theory should apply (Lumley and Panofsky, 1964). An important aspect of this is that the dimensionless flux profiles are universal functions, , of the stability coordinate ζ = d/L, where L is the Monin-Obukhov length scale,
k is the von Karman constant, B0 is the surface buoyancy flux into the ocean, and (u*)2 is the kinematic surface wind stress. This leads to logarithmic mean-property profiles near the boundary; the associated fluxes are about 20 percent of their surface values at and approach their surface values linearly as d→0. Since there is no turbulent transfer either across the surface, or at d = h if the interior mixing is zero, the diffusivities must satisfy Kx(0) = Kx(h) = 0. The universal shape is assumed to be cubic in d, which has had rough empirical confirmation; thus,
The turbulent velocity scale wx is formulated for consistency with similarity theory as a simple combination of the surface forcing velocities u* and w* (which is equal to (−B0h)1/3 if B0 is less than 0).
The counter-gradient term γx is non-zero only for the buoyancy field and for passive scalars in unstable forcing conditions (i.e., B0 < 0); thus, we restrict attention to x = s (for scalar). As suggested by Deardorff (1972), it has been successfully parameterized as
with Cs = 6.2.
The thickness of the oceanic planetary boundary layer h is determined as the smallest value of d at which a bulk Richardson number Rib across the boundary layer achieves a certain critical value Ric. Thus it primarily depends on the buoyancy profile b(d) and the velocity profile v(d), but there is also a (usually weak) dependence on the turbulent intensity:
The turbulent velocity vt combines with the mean velocity difference in the denominator of (14); its primary role is to establish the correct entrainment flux at the interior edge of the PBL when the mean velocity is weak or absent (i.e., primarily in the case of free convection). There is strong empirical evidence that this entrainment buoyancy flux is a fixed fraction, βT, of the surface buoyancy flux B0 = − , viz., βT = 0.20 (see, e.g., Tennekes, 1973; this result is sometimes known as Turner's rule). With the parameterization forms above, it is achieved for
where ct = 5.8 is a constant that arises from the empirically determined stability function . Thus, from (16) and (17), h is larger, and hence the entrainment layer is thicker, where N(h) is smaller.
We have used this model, together with a new parameterization of internal vertical mixing described below that joins smoothly onto the PBL mixing at d = h, to calculate a variety of solutions, on time scales from hours to decades. Among these is one coupled to a one-dimensional radiative-convective atmosphere for the Ocean Weather Station Papa site in the North Pacific Ocean (Figures 9 through 11). The quality of its annual-cycle simulation is very good by the standards of previous upper-ocean PBL models (Martin, 1985), and we have found that an important part of this success is the inclusion of a diurnal cycle in insolation and an idealized, periodic storm cycle in the surface winds. To assess the consequences for climate modeling of accurately representing the oceanic annual cycle, we can compare Figure 11 with another coupled solution whose "ocean" is no more than a constant-depth mixed layer; it shows differences in tropospheric temperatures of up to 5 K in the annual cycle.
We are now in the process of extending and calibrating this new parameterization using large-eddy simulation solutions for the PBL under a variety of stress and buoyancy forcing conditions (e.g., Wyngaard and Moeng, 1993; McWilliams et al., 1993). We are particularly concerned with representing the asymmetry of scalar diffusivities forced by either surface or entrainment fluxes under convective conditions (Wyngaard and Brost, 1984). This effect is absent from the present parameterization forms, but it may be quite important for the diverse oceanic combinations of