\documentstyle [12pt]{article}

\oddsidemargin 0.0in     \evensidemargin 0.0in
\topmargin -0.5in        \headsep .5in
\textheight 8.8in        \textwidth 6.5in
\parskip 7pt

\def\baselinestretch{1.5}

\begin{document}

\vspace*{3cm}

\begin{center}
\bf
Mixed-Layer/Thermocline Interaction in a Three-Dimensional Isopycnic 
Coordinate Model \\
\vspace{3cm}
\rm
by \\
\vspace{3cm}
Rainer Bleck\footnotemark[1], Howard P. Hanson\footnotemark[2], 
Dingming Hu\footnotemark[1], Eric B. Kraus\footnotemark[2] \\
\vspace{2cm}
Revised Version -- February 1989
\end{center}

\footnote[1]{Rosenstiel School of Marine and Atmospheric Sciences, 
University of Miami, Miami, FL 33149}
\footnote[2]{Cooperative Institute for Research in Environmental Sciences, 
University of Colorado/NOAA, Boulder, CO 80309-0449}

\newpage

\begin{center}

ABSTRACT

\end{center}

The annual buildup and obliteration of the seasonal thermocline and the 
associated ventilation of the permanent thermocline in a wind- and thermally- 
driven ocean basin are simulated numerically. The model developed for this 
purpose is a combination of a single-layer model of the oceanic mixed layer, 
based on a simple closure of the turbulence kinetic energy equation, and a 
three-dimensional isopycnic coordinate model of the stratified oceanic 
interior. The joint model, set in a rectangular ocean basin, is forced by 
annually varying wind stress and radiative plus turbulent heat fluxes 
approximating zonally averaged conditions over the North Atlantic. Special 
emphasis is placed on the description of the mixed layer detrainment process, 
which requires distributing mixed-layer water of continuously variable 
density among constant-density interior layers. The truncation errors 
associated with this process are found to be numerically tolerable.

The quasi-Lagrangian character of the model's vertical coordinate permits 
easy tracking of water masses that are left behind during the annual retreat 
of the mixed layer to form the seasonal thermocline. Likewise, the subduction 
of ventilated water into the permanent thermocline by the horizontal gyre 
motion is explicitly simulated. 

While a comparison of simulated mixed-layer characteristics with actual 
observations is problematic due to the idealized basin configuration, the 
model appears to be reasonably successful in duplicating the seasonal 
cycle of the zonally averaged conditions over the North Atlantic.

\newpage

1. {\bf Introduction}

The response of the ocean to changes in atmospheric conditions, both 
man-made and natural, is one of the central topics of climate research 
today. Processes that determine the thermal structure of the thermocline, 
its motion field, and its buffering capacity are among the topics being 
studied. Many of these investigations have been inspired by the ventilated 
thermocline model of Luyten {\em et al.} (1983), who succeeded in deriving 
analytic solutions for the flow in interior isopycnic layers forced by 
density and potential vorticity boundary conditions at the bottom of the 
Ekman layer in regions of downward Ekman pumping. Many features of their 
analytic solutions (for example, the existence of unventilated ``shadow" 
and ``pool" zones) are supported by observations and have been reproduced 
in numerical simulations by Cox and Bryan (1984), Cox (1985), and Huang 
and Bryan (1987).

One element common to the above-mentioned thermocline models is the use of 
time-independent external forcing functions. Woods (1985) has emphasized 
the need to incorporate the annual growth and retreat of the oceanic mixed 
layer into these models, arguing that it is this process that determines the 
potential vorticity distribution in the ventilated thermocline. We report 
here on an effort to simulate the annual mixed layer entrainment/detrainment 
cycle in a three-dimensional numerical circulation model that combines 
isopycnic coordinate representation of the interior adiabatic flow with a 
simple mixed-layer model driven by seasonally varying buoyancy and wind 
forces. While the emphasis in our modeling work has been on the mass exchange 
between mixed layer and oceanic interior and on the gyre-scale migration of 
subducted water masses, we intend to demonstrate that the resulting model 
does a reasonable job of simulating the annual cycle of sea surface temperature 
and upper-ocean heat content.

It is our intention to equip an existing isopycnic model of the North and 
Equatorial Atlantic with a similar type of mixed layer parameterization; 
however, in this paper we will discuss results obtained in an idealized, 
rectangular, flat-bottom ocean basin forced by the annual cycle of zonally 
averaged wind and buoyancy fluxes approximating the climatic conditions 
over the Atlantic Ocean between $ 10^{0} $ and $ 60^{0} $ N.

The model presented here offers some advantages over conventional Cartesian 
coordinate models with regard to depicting subduction and thermocline 
ventilation processes:

(1) At every horizontal grid point the depth of the mixed layer is 
unambiguously defined at all times. This allows us to draw a clear 
distinction between mixed layer and interior water.

(2) The amount of water discharged from the mixed layer into a particular 
isopycnic layer is also explicitly computed. This water can be ``tagged" 
numerically and can be followed on its path around the wind-forced gyres. 
In particular, it is possible to keep track of that part of the detrained 
water that escapes the following year's re-entrainment into the mixed 
layer, i.e., the water that migrates from the seasonal into the permanent 
thermocline.

Isopycnic coordinate models offer a solution to some, though certainly not 
all problems encountered in Cartesian coordinate models. Since vertical 
resolution is of necessity finite, truncation errors will, in one way 
or another, limit the realism of the numerical solution. Judging from 
shortcomings of other types of mixed-layer models, these errors should 
come to light primarily during mixed-layer detrainment. And indeed they 
do. The specific problem posed by the isopycnic framework is that there 
are only a finite number of isopycnic layers available to accept the 
detrained water (which, as long as it resides in the mixed layer, is 
allowed to have a density varying in space and time). Some errors in the 
timing and the rate of detrainment, and consequently in the potential 
vorticity of the detrained water masses, are therefore inevitable if the 
model is run in very low vertical resolution mode. How significant these 
errors are compared to the truncation errors encountered in Cartesian 
coordinate models using the same number of grid points in the vertical is 
difficult to assess. There is reason to be believe, though, that a Cartesian 
model with a comparable number of grid points in the vertical would be 
unable to simulate the annual cycle of mass exchange between the thermocline 
and the mixed layer with the same degree of realism as the present model.

Since we are exploring new modeling technology we consider it wise to 
conduct our first experiments in low resolution mode, both horizontally and 
vertically. In particular, the present model is not intended to resolve 
scales of motion typical for barotropic and baroclinic instability in the 
ocean. Vertical resolution is also marginal at best. Due to the 
uncertainties surrounding the magnitude and numerical representation of 
diapycnic mixing processes we have decided to omit those from the model 
for the time being. The results of this omission are as expected and are 
considered tolerable: at high latitudes the model continues to produce 
relatively dense water, which accumulates in the basin at the rate of a 
few meters per year. To reach a steady state in the water mass distribution, 
the model would have to be run long enough for the ``Deep Water" to reach 
the surface at low latitudes where it could then be heated and recycled 
through the mixed layer. We have not attempted to carry the model 
integration to this unrealistic final stage. 

The paper is structured as follows: After providing a model description in 
Sections 2 and 3 we will discuss in Section 4 our method of selecting the 
wind stress and thermal forcing functions. Basic features of the resulting 
circulation will be presented in Section 5. In Section 6 we will discuss 
annual mixed layer variations in greater detail. The topic of Section 7 
will be the migration of ventilated water in the permanent thermocline. 
Section 8 provides a summary and suggestions for future experiments.

\vspace{2cm}

2. {\bf Model Layout, Basic Equations}
\nopagebreak

The dynamic equations of the model are solved on a grid consisting of 
32x32 horizontal grid points and 6 layers in the vertical. North-south 
grid lines correspond to meridians spaced $ 2^{0} $ apart, yielding a total 
basin width of $ 64^{0} $ longitude. East-west grid lines run at constant 
latitude. The grid points form a regular mesh on a Mercator projection, 
i.e., the north-south spacing of grid points on the sphere changes with 
the spacing of the meridians. The computational domain ranges in latitude 
from $ 11.9^{0} $ to $ 59.3^{0} $. The mesh size varies from 218 km in the 
south to 114 km in the north. The basin depth is 5 km.

Of the six layers in the vertical, the uppermost is occupied by the 
thermodynamically active mixed layer while layers 2 ... 6 are layers of 
constant specific volume $ \alpha = \rho^{-1} $. Their assigned coordinate 
values, expressed here in non-dimensional units of $ (1 - \alpha/\alpha_{0}) 
\times 10^{3} $, are $ \sigma = $ 26.1, 26.5, 26.9, 27.3, 27.7 where 
$ \alpha_{0} = 10^{-3} m^{3}/kg. $ (We use the symbol $ \sigma $ since 
for practical purposes these coordinate values are identical 
to conventional ``$ \sigma $ units", $ (\rho/\rho_{0} - 1) \times 10^{3} $ 
with $ \rho_{0} = \alpha_{0}^{-1}. $) All isopycnic layers, even those that 
contain mass only in limited geographical regions, are carried by the model 
at every horizontal grid location. This arrangement does not undermine static 
stability as long as layers whose $ \sigma $ value falls below the spatially 
and temporally variable $ \sigma $ value in the mixed layer have zero mass. 

The prognostic equations being solved in the isopycnic interior are 
identical to those used by Bleck and Boudra (1986) in their ``pure" 
isopycnic model. The set, repeated here for convenience, consists of the 
layer-integrated horizontal momentum equation
\[
\frac{\partial \vec{v}}{\partial t} + \nabla \frac{\vec{v}^{2}}{2} 
+ (\zeta + f)\vec{k} \times \vec{v} + \nabla M = \frac{1}{\Delta p} 
\nabla \cdot (\nu \Delta p \: \nabla \vec{v})
\]
and the layer-integrated continuity equation
\[
\frac{\partial}{\partial t} \Delta p + \nabla \cdot (\vec{v} \Delta p) = 0.
\]
Here, $ \vec{v} = (u,v) $ is the horizontal velocity vector, 
$ \zeta = (\partial v / \partial x)_{\alpha}
        - (\partial u / \partial y)_{\alpha} $ 
is the  relative vorticity, $ M \equiv  gz + p \alpha $ is the Bernoulli 
function or Montgomery potential, $ gz $ is the geopotential, $ \Delta p $ 
is the pressure thickness  of the isopycnic layer, $ f $ is the Coriolis 
parameter, $ \vec{k} $ is the vertically pointing unit vector, and $ \nu $ 
is a variable eddy viscosity coefficient proportional to the large-scale 
deformation (Smagorinsky, 1963). The layer-depth weighting in the viscosity 
term is necessary to conserve overall momentum. All horizontal derivatives 
are evaluated at $ \alpha = constant $ and the local time derivative applies 
to a fixed point in $ x,y,\alpha $ space. Distances in $ x $ and $ y $ 
direction are measured in the projection onto a horizontal plane, implying 
that the velocity components $ (u,v) \equiv (dx/dt,dy/dt) $ are the 
"horizontal" components of motion in the literal sense. As discussed by 
Bleck (1978), these conventions prevent the generation of so-called metric 
terms related to the slope of the $ \alpha $ surface in deriving the above 
equations. 

Under adiabatic flow conditions, which we assume to prevail in the isopycnic 
interior, the thermodynamic equation is satisfied by setting the vertical 
velocity component $ d \alpha /dt $ to zero, i.e., by treating the isopycnic 
layer interfaces as material surfaces. Note that this simplification has 
already been incorporated into the continuity equation shown above.

The set of equations to be solved in the isopycnic interior is completed 
by the hydrostatic equation
\[
\frac{\partial M}{\partial \alpha} = p.
\]

The equations solved in the mixed layer differ from the above set in three 
respects:

(a) The horizontal pressure force term in the momentum equation, which in 
isopycnic layers equals $ - \nabla M $, must be expanded into a two-term 
expression $ - (\nabla M  - p \nabla \alpha) $ to account for the fact that 
the mixed layer is non-isopycnic.

(b) A wind stress term must be added to the right-hand-side of the momentum 
equation. If we stipulate that the direct effect of the the wind forcing 
does not extend past the bottom of the mixed layer, which implies that the 
wind-induced stress profile $ \tau(z) $ reaches zero somewhere between 
the sea surface and the mixed layer depth $ \Delta p $, then the stress 
term in the layer-integrated momentum equation assumes the form 
$ (\Delta p)^{-1} g \vec{\tau}(0) $. (In this context it is worth noting 
that there are no interfacial stresses among coordinate layers. The 
interior layers are set in motion exclusively by pressure torques and, to a 
lesser extent, by bulk momentum transfer during mixed layer detrainment.)

(c) In case of material layer interfaces, the change in layer thickness 
$ \partial (\Delta p) / \partial t $ brought about by con- or diverging 
mass fluxes within the layer is synonymous with the vertical velocity 
difference of water particles at the upper and lower interface. Since the 
bottom of the mixed layer is not a material surface, the material vertical 
motion of water parcels at the bottom of the mixed layer given by $ \partial 
(\Delta p) / \partial t $ must be split into two components (Bleck, 1978): 
the actual (nonmaterial) vertical displacement of the coordinate surface 
marking the bottom of the mixed-layer and the movement of water $ relative $ 
to this surface. The latter quantity is the actual entrainment/detrainment 
rate. The discussion of its diagnosis will occupy a major part of Section 3. 

Since the mixed layer is non-isopycnic and subject to diabatic forcing, 
conservation equations for thermodynamic variables are needed there. In the 
present model version, salinity is held constant at 34.5 mil. This reduces 
the required set of thermodynamic equations to a single transport equation 
for mixed layer temperature: 
\[
\frac{\partial T}{\partial t} + \vec{v} \cdot \nabla T = g \frac{H+R}
{\Delta p c_{w}}.
\]
Here, $ H $ and $ R $ are the turbulent and radiative heat fluxes through 
the sea surface and $ c_{w} $ is the specific heat of sea water. 
(Entrainment fluxes are not included here and will be discussed later.)
For the equation of state relating temperature $ T $ (in $ ^{0}C $) and 
salinity $ S $ (in mil) to $ \sigma $ we use a power series expansion 
quadratic in $ T $ and linear in $ S $:
\[
\sigma = c_{0} - c_{4} \left( \frac{c_{2}}{2} (T^{2}-25) + c_{5}(T-5) 
- (S-35) \right)
\]
where $ c_{0} = 27.67547 $, $ c_{1} = 0.07 $, $ c_{2} = 0.013 $, 
$ c_{3} = 0.004 $, $ c_{4} = 0.0008 $, and $ c_{5} = c_{1} + c_{3}(S-35) $. 
This equation of state, a simplified version of the equation used by 
McDougall (1987), has the advantage of giving a closed expression for 
temperature as a function of density and salinity:
\[
T = \frac{1}{c_{2}} \left[ -c_{5} + \sqrt{ (c_{5} + 5c_{2})^{2} 
+ 2 c_{2} \left( \frac{c_{0}-\sigma}{c_{4}} + (S-35) \right)} \right].
\]
With this formula we can compute the temperatures corresponding to the 
five $ \sigma $ values assigned to the isopycnic coordinate layers. For 
$ S = 34.5 $ we obtain the values 13.0, 10.7, 8.1, 4.8, $ -0.4^{0} C $.

The model equations are solved using standard finite-difference techniques 
which we will not reproduce in full detail here. The following items seem 
worth mentioning, however.

The equations are solved on a ``B" grid, which carries the velocity components 
$ u,v $ at the four corners of each grid box and mass field variables 
($ p, T, \sigma, M $) in the center. The horizontal mass-flux components 
$ u \Delta p, v \Delta p $ are defined halfway between mass points in $ x $ 
and $ y $ direction respectively to reduce the truncation error in 
computing $ \nabla \cdot (\vec{v} \Delta p) $. The mass-flux curl, which 
is needed to filter out the barotropic mass flux divergence (see below) 
is then most appropriately defined at $ u,v $ points. The actual side wall 
of the ocean basin coincides with $ u,v $ points. Both velocity components, 
as well as the streamfunction for the mass flux, are set to zero there 
(no-slip boundary condition).

Variables are also staggered in the vertical: $ p $ and $ gz $ are referred 
to as ``level" variables and are attached to layer interfaces. The rest, 
i.e., $ u, v, T, \sigma, M $, are ``layer" variables which are 
height-independent within each coordinate layer. The layer thickness 
$ \Delta p $ is also considered a layer variable.

The numerical scheme for solving the continuity equation in partially 
"collapsed" isopycnic layers (i.e., in the presence of zero $ \Delta p $ 
values) is discussed in detail in Bleck and Boudra (1986). 

The ``B" grid is known for its tendency to produce checkerboard noise in the 
mass field. To suppress this noise we smooth the layer interface pressure 
field after each time step. This is done by replacing the $ p $ value on 
each interface at each horizontal grid point by $ p + \nabla \cdot (\nu 
\nabla p) \Delta t $ where $ \nu = 1.5 cm\:s^{-1} \times (mesh \; size) $. 
The ``diffusive pressure flux" $ \nu \nabla p $ across lateral boundaries 
is set to zero. One side effect of this operation is that interface slopes 
are reduced: the mass field becomes less baroclinic. The fluid reacts by 
decreasing the vertical velocity shear. Pressure smoothing therefore entails 
vertical momentum mixing. Such transfer of momentum from baroclinic to 
barotropic mode is considered beneficial because it duplicates some effects 
of baroclinic instability.

Caution must be exercised in formulating the horizontal pressure force 
term in the mixed layer, $ \nabla M - p \nabla \alpha $ = $ \nabla gz 
+ \alpha \nabla p $, to assure that pressure torques are transmitted 
correctly across the interface between mixed layer and isopycnic interior. 
If those torques, whose magnitude strongly depends on interface tilt, are 
not correctly passed from layer to layer, the finite difference 
approximation to the curl of the vertically integrated pressure force 
will not reduce to a simple bottom torque term and the basic Sverdrup 
balance will be disturbed. Especially during winter, when the mixed layer 
bottom exhibits a strong meridional tilt, serious errors in gyre strength 
can result. (The fluctuations in subpolar gyre strength shown in Fig. 5 of 
Kraus {\em et al.}, 1988, are due to this problem). 

The finite difference expressions 
\[
\frac{\partial gz}{\partial x} + \alpha \frac{\partial p}{\partial x} = 
\overline{\delta_{x} \overline{gz}^{s} + \frac{\overline{\alpha \Delta p}
^{x}} {\overline{\Delta p}^{x}} \delta_{x} \overline{p}^{s}}^{y}
\]\[
\frac{\partial gz}{\partial y} + \alpha \frac{\partial p}{\partial y} = 
\overline{\delta_{y} \overline{gz}^{s} + \frac{\overline{\alpha \Delta p}
^{y}} {\overline{\Delta p}^{y}} \delta_{y} \overline{p}^{s}}^{x}
\]
have very nearly the desired property. Here, overbar-x, overbar-y, and 
overbar-s define two-point averaging operations in $ x $, $ y $, and 
vertical direction respectively, while $ \delta_{x} $ and $ \delta_{y} $ 
denote two-point, centered finite differences in the two horizontal 
directions. The salient element of the above expressions is the 
horizontal layer-depth weighting of $ \alpha $.

Kinetic energy is dissipated at the lateral boundaries through the use of 
no-slip boundary conditions ($ u,v = 0 $). An additional sink of kinetic 
energy is provided by a bottom drag of the form $ - \rho c_{D} |\vec{v}| 
\: \vec{v} $ where $ \vec{v} $ is the velocity in the 6th layer and $ c_{D} 
= 0.003 $.

In order to extend the maximum permissible time step, barotropic gravity 
waves are removed from the numerical solution by correcting the $ u,v $ 
field after each time step to eliminate the barotropic mass flux 
divergence. While this filtering device is commonly referred to as a rigid 
lid, in the present non-Boussinesq model it actually has the effect of 
keeping the bottom pressure time-independent. In connection with a leap 
frog time integration scheme, this rigid-lid equivalent allows us to extend 
the time step to 1/15 of a day.

The earth's sphericity is taken into account by incorporating a variable 
map scale factor into the finite-difference equations. Additional metric 
terms in the equations are avoided by formulating vorticity and divergence 
as line integrals around individual grid boxes.

The equations presented so far make no reference whatsoever to the exchange 
of mass between the mixed layer and the interior layers, ostensibly one 
of the central themes of our modeling effort. In building these exchange 
processes into the model framework we have found it expedient not to 
represent them as additional source terms in the above equations. Rather, 
the mass exchange is being treated numerically as an instantaneous adjustment 
to the vertical mass and momentum distribution in each grid column at the 
end of each prognostic time step. Only in this way are we able to accomodate 
the diversity of layer thickness configurations, including the presence 
of one or more ``massless" isopycnic layers beneath the mixed layer, without 
being overwhelmed by mass, momentum, and thermal energy conservation 
requirements, not to mention non-negative layer thickness constraints. This 
is the topic of the next section.

\vspace{2cm}

3. {\bf The Mixed Layer}
\nopagebreak

3a. {\em Formulation of the Mixed Layer Energy Balance}
\nopagebreak

As outlined in the original theoretical development of the ``Integral 
Model'' by Kraus and Turner (1967), the distinguishing feature of the 
mixed layer, in contrast to the oceanic interior, is that motion in it 
is sufficiently turbulent to achieve complete homogenization. Turbulence 
is maintained against viscous dissipation by wind stress and cooling 
from above (upward buoyancy flux). Excess turbulence kinetic ernergy is 
postulated to be spent instantaneously on extending the range of the 
turbulent region, i.e., on generating potential energy by mixing heavier 
water upward from the stratified oceanic interior.

In the simplest implementation of the Integral Model (e.g., Niiler and 
Kraus, 1977), turbulence kinetic energy (TKE) input into the mixed layer 
through wind stress is parameterized in terms of $ u^{*3} $, the third 
power of the friction velocity. A fixed percentage of this energy is 
assumed to be lost to viscous dissipation. The remainder, i.e., the net 
gain in TKE due to wind forcing, is traditionally expressed as $ m u^{*3} 
$ where $ m $ is an empirically derived parameter. (We use $ m = 1.25 $.)

Depending on its sign, the vertical buoyancy flux in the mixed layer acts 
either as a sink or source of TKE. In the latter case (which is realized if 
the buoyancy flux is directed upward, i.e. if the ocean loses heat to the 
atmosphere), a fixed percentage of TKE generation is again assumed to be 
lost to dissipation. The residual buoyancy flux available for TKE 
generation, averaged over the mixed layer, is traditionally expressed 
as $ n B_{0}/2 $ where $ n < 1 $ is another empirical parameter and 
$ B_{0} $ is the positive (upward) buoyancy flux at the sea surface. (We 
use $ n = 0.4 $.) The factor 2 in the denominator reflects the assumption 
that, in circumstances of maximum dissipation, the buoyancy flux varies 
linearly between zero at the bottom of the mixed layer and $ B_{0} $ at 
the top. Since buoyancy is the product of gravity and relative density 
variation, $ B_{0} $ is dimensioned $ m^{2}/s^{3} $. Note that it is 
the buoyancy flux itself, rather than its divergence, that appears as a 
forcing term in the TKE equation. To express it in the form of an energy 
flux through the sea surface, in analogy to the wind stress term the 
buoyancy flux must therefore be integrated over the mixed layer depth 
$ h $. The resultant energy generation term in case of positive (upward) 
buoyancy flux is $ n h B_{0} / 2 $.

When the ocean is heated from above, the buoyancy flux is directed downward 
and acts to suppress TKE in the mixed layer. By again assuming that the 
buoyancy flux decreases linearly to zero over the depth of the mixed layer, 
$ B_{0} < 0 $ implies a negative (outgoing) TKE flux through the surface of 
the form $ h B_{0}/2 $. In summary, if $ B_{0} $ is positive, it contributes 
positively to TKE generation at the rate $ n h B_{0}/2 $, while a negative 
$ B_{0} $ contributes negatively to TKE generation at the rate $ h B_{0}/2 $.

As mentioned earlier, any net TKE increase is assumed to be spent on mixing 
up heavier water from below the interface. If $ b_{mix} $ denotes the 
buoyancy in a mixed layer of depth $ h $, and $ b_{sub} $ denotes the 
buoyancy immediately below the mixed layer, then the energy needed to mix 
a slab of thickness $ dh $ into the mixed layer is given by $ (b_{mix} 
- b_{sub}) h dh/2 $. The rate at which this process soaks up kinetic 
energy is
\begin{equation}
(b_{mix}- b_{sub}) w_{e} \frac{h}{2} \label{eq:1}
\end{equation}
where $ w_{e} = dh/dt $ ($ w_{e} = dh/dt - w $ in case of nonzero mean 
vertical motion $ w $) is the entrainment velocity. Equating this to the 
net TKE generation rate yields
\[
m u^{*3} + \left( \frac{B_{0} - |B_{0}|}{2} + 
                n \frac{B_{0} + |B_{0}|}{2} \right) \frac{h}{2} + 
(b_{mix}- b_{sub}) w_{e} \frac{h}{2} = 0
\]
where the two possibilities of upward and downward buoyancy flux have been 
combined into one expression. It follows that the entrainment rate can be 
obtained from the expression
\begin{equation}
w_{e} = \frac{(2/h) m u^{*3} + [( B_{0} - |B_{0}|) + 
                              n ( B_{0} + |B_{0}|)] / 2}
{b_{mix}- b_{sub}}. \label{eq:2}
\end{equation}
This formula obviously loses validity if it yields $ w_{e} < 0 $, which 
happens if the stabilizing effect of a negative buoyancy flux outweighs 
the stress-related positive TKE source term. In such a case of negative 
TKE generation Kraus and Turner (1967) assert that turbulence cannot be 
maintained over the depth range $ h $ and that $ h $ should be adjusted 
until (\ref{eq:2}) yields equilibrium, i.e., $ w_{e} = 0 $. The mixed layer 
depth where this occurs is
\[
h = \frac{2 m u^{*3}}{ -B_{0}}
\]
which, for $ m = 1.25 $, agrees with the conventional definition of the 
Monin-Obukhov length, $ L = u^{*3} / (- \kappa B_{0}) $. This value for 
$ m $ turns out to agree well with experimental results (Niiler and Kraus, 
1977).

3b. {\em Numerical Implementation of Mixed-Layer Entrainment/Detrainment}
\nopagebreak

The present model consists of a finite number of isopycnic layers marked 
by discrete buoyancy values  $ b_{k} = - g \sigma_{k} $, capped by a 
mixed layer in which $ b $ is a function of $ x,y,t $. Not all isopycic 
layers in a given grid column are required to contain mass; in fact, 
layers whose assigned buoyancies $ b_{k} $ exceed the mixed layer value 
$ b(x,y,t) $ must be in a ``massless" state to prevent static instability. 
In the course of the annual heating/cooling cycle those upper isopycnic 
layers may periodically be ``inflated" by water detrained from the mixed 
layer and thus act as receptacle for water forming the seasonal thermocline.

While a deepening mixed layer is relatively easy to model numerically, the 
detrainment process poses a far greater challenge, no matter what type of 
mixed layer model one uses. As alluded to in the Introduction, the main 
difficulty in the present model is as follows. The buoyancy $ b $ of the 
detrained column generally does not match the buoyancy of an isopycnic 
layer but falls between the buoyancies of two coordinate layers, say, 
$ b_{k} $ and $ b_{k-1} $. Before part of the mixed layer can be 
"backpumped" into existing isopycnic layers, thermal energy must therefore 
be redistributed within the mixed layer with the goal of either

(a) reducing the buoyancy in the column to be detrained to $ b_{k} $; or

(b) reducing the buoyancy in one part of the column to be detrained to 
$ b_{k} $ and raising it in the remaining part to $ b_{k-1} $.

In both cases heat must be extracted from the detrained column and 
deposited in the upper, remaining part of the mixed layer. In case (a) this 
is done for the purpose of conserving overall thermal energy, while in case 
(b) the heat transfer is needed to maintain a stable stratification in the 
column. 

In a model relying on bulk parameterization of surface heat fluxes, a 
numerically induced increase in surface temperature during the detrainment 
phase may have undesirable consequences. We will discuss this further in the 
section dealing with the detrainment algorithm.

3c. {\em The Entrainment Algorithm}
\nopagebreak

Integrating (\ref{eq:2}) over a time step $ \Delta t $ yields a formula for 
the depth $ \Delta h $ of the water column entrained into the mixed layer 
during $ \Delta t $:
\begin{equation}
\Delta h = \frac{E \Delta t}{b_{mix}- b_{sub}} \label{eq:3}
\end{equation}
where
\[
E \Delta t = \left[ \frac{2}{h} m u^{*3} +
   \frac{B_{0} - |B_{0}|}{2} +
 n \frac{B_{0} + |B_{0}|}{2} \right] \Delta t.
\]
Note that (\ref{eq:3}) cannot be incorporated into a model in its present 
form because the denominator is not bounded away from zero. Furthermore, 
(\ref{eq:3}) is built upon the assumption that water of buoyancy $ b_{sub} $ 
actually exists throughout the depth range $ \Delta h $. In order to avoid 
possible inconsistencies (or outright numerical overflow) in diagnosing 
$ \Delta h $ from (\ref{eq:3}) it is advisable to replace the expression 
(\ref{eq:1}) by a more general one that takes into account the actual 
isopycnic layer structure underneath the mixed layer and yields an exact 
expression for the number of layers the mixed layer will entrain before the 
finite energy supply $ E \Delta t $ is exhausted. Such an alternate 
expression is obtained as follows.

Let $ b_{1} $ and $ z_{1} $ denote the buoyancy and the depth of the mixed 
layer at a particular horizontal location and let $ (z_{2}, z_{3}, ..., 
z_{n}) $ denote the lower interfaces of isopycnic layers $ (b_{2}, b_{3}, 
..., b_{n}) $, any number of which may be massless in the grid column 
considered. The potential energy of the above ensemble of layers prior 
to mixing is then given by
\begin{equation}
\int_{z_{0}}^{z_{n}} b z dz = \frac{1}{2} 
\sum_{k=1}^{n} b_{k} ( z_{k}^{2} - z_{k-1}^{2} ), \label{eq:4}
\end{equation}
whereas the potential energy in the same column after vertical mixing is
\begin{equation}
b_{mix} \int_{z_{0}}^{z_{n}} z dz = b_{mix} \frac{z_{n}^{2}}{2} \label{eq:5}
\end{equation}
with
\begin{equation}
b_{mix}= \frac{1}{z_{n}} \sum_{k=1}^{n} b_{k} ( z_{k} - z_{k-1} ) 
\label{eq:6}
\end{equation}
and $ z_{0} = 0 $. Note that averaging the individual buoyancy values to 
obtain $ b_{mix} $ as shown in (\ref{eq:6}) is only permitted if 
nonlinearities in the equation of state for sea water can be ignored. These 
nonlinearities are such that (\ref{eq:6}) slightly underestimates the value 
of $ b_{mix} $. Linearization of the equation of state in evaluating 
(\ref{eq:6}) seems justified in light of the uncertainties involving the 
constants $ m $ and $ n $ in the TKE generation term.

Combining (\ref{eq:5}) and (\ref{eq:6}) yields the following expression for 
the potential energy of the column after mixing:
\begin{equation}
b_{mix} \int_{z_{0}}^{z_{n}} z dz = 
\sum_{k=1}^{n} b_{k} ( z_{k} - z_{k-1} ) \frac{z_{n}}{2}. \label{eq:7}
\end{equation}
If the difference between (\ref{eq:7}) and (\ref{eq:4}) were equal to 
$ E \Delta t $, the interface depth $ z_{n} $ would be the sought-after 
entrainment depth $ h $. Since this is not likely to happen, we must 
substitute $ h $ for the interface depth $ z_{n} $ in (\ref{eq:7}) and 
(\ref{eq:4}) and treat $ h $ as the unknown in the problem. 

Equating the difference (\ref{eq:7})--(\ref{eq:4}) to $ E \Delta t $ and 
solving for $ h $ gives
\begin{equation}
h = \frac{2 E \Delta t + G(n-1) - b_{n} z_{n-1}^{2}}
{F(n-1) - b_{n} z_{n-1}} \label{eq:8}
\end{equation}
where
\[
F(n) = \sum_{k=1}^{n} b_{k} ( z_{k} - z_{k-1} )
\]
and
\[
G(n) = \sum_{k=1}^{n} b_{k} ( z_{k}^{2} - z_{k-1}^{2} ).
\]
Before $ h $ can be computed, $ n $ must be known. In practice, both tasks 
are solved simultaneously by computing ``trial" values of $ h $ successively 
for $ n = 2, 3, ... $. This is continued as long as (\ref{eq:8}) yields an 
entrainment depth greater than the upper interface of the $ n $-th layer, 
$ z_{n-1} $. In other words, the correct value of $ n $ and the correct 
entrainment depth become known as soon as the process of incrementing 
$ n $ has gone one step too far. 

Note that the above algorithm, which comes to completion after a finite 
number of steps and is readily vectorized, properly accounts for massless 
isopycnic layers sandwiched between the mixed layer and the ocean interior.

3d. {\em The Detrainment Algorithm}
\nopagebreak

The entrainment process described in the previous section is an 
"instantaneous" rearrangement of the vertical density profile in response 
to TKE input during the preceding time interval $ \Delta t $. Likewise, the 
detrainment process described here is an instantaneous transfer of mass 
and thermal energy among various model layers that takes place at the 
end of the time interval $ \Delta t $ when the TKE input $ E \Delta t $ 
is negative.

Consider the situation at the end of time interval $ \Delta t $ when buoyancy 
and mass in the mixed layer have been transported horizontally and buoyancy 
has been modified by thermal energy fluxes through the surface, but before 
detrainment has taken place. As before, use subscript $ 1 $ to identify 
mixed layer buoyancy and depth, and denote isopycnic layer densities and 
lower interface depths by subscripts $ 2,3,...$. The corresponding values 
after detrainment will be identified by primes. Let $ b_{k} $ and $ b_{k-1} $ 
be the two layer buoyancies that bracket $ b_{1} $, and let $ L $ be the 
Monin-Obukhov length computed from the local stress and buoyancy fluxes. 
$ L $ must be smaller than $ z_{1} $ to set the stage for detrainment.

The simplest case to consider is the one in which the total column of 
detrained water, $ z_{1} - L $, is cooled until it matches the buoyancy 
$ b_{k} $. The thermal energy to be transferred to the water column above 
depth $ L $ will then raise the buoyancy there to
\begin{equation}
b_{1}' = b_{1} + ( b_{1} - b_{k} ) \frac{z_{1} - L}{L}. \label{eq:9}
\end{equation}
(As in the entrainment case we ignore nonlinearities in the equation of 
state for the time being, that is, we redistribute buoyancy as if it were 
thermal energy. More will be said about this subject later.) 

The buoyancy transfer described in (\ref{eq:9}) concludes the detrainment 
process if $ both $ of the following conditions are met:

Condition A: The value $ b_{1}' $ given by (\ref{eq:9}) falls below a 
physically realistic threshold value $ b_{max} $, a variable meant to 
safeguard against errors in surface heat flux calculations resulting from 
an ``overheated" mixed layer. The logical choice for $ b_{max} $ is the 
buoyancy the mixed layer would have obtained during the preceding time 
interval $ \Delta t $ had the incoming heat been distributed only over 
the ``new" mixed-layer depth $ L $ rather than the ``old" depth $ z_{1} $. 
This leads to the following expression for $ b_{max} $:
\[
b_{max} = b_{1} + (-B_{0}) \Delta t  ( \frac{1}{L} - \frac{1}{z_{1}} ).
\]

Condition B: The value $ b_{1}' $ falls below $ b_{k-1} $.

If condition B is violated, the necessity for detraining mixed layer water 
exclusively into isopycnic layer $ k $ no longer exists. Water can and 
should be divided up among layers $ k $ and $ k-1 $ in this case. However, 
before describing this variant of the detrainment scheme the implications 
of condition A need to be discussed.

If $ b_{1}' $ exceeds $ b_{max} $, part of the water set aside for 
detrainment into layer $ k $ needs to be retained to lower the mixed layer 
buoyancy to $ b_{max} $. The net effect of this adjustment is that the 
mixed layer, instead of receding to $ L $, recedes from $ z_{1} $ to the 
depth
\begin{equation}
z_{1}' = z_{1} \frac{b_{1}-b_{k}}{b_{max} - b_{k}}. \label{eq:10}
\end{equation}

Having found a way to prevent mixed layer buoyancy from exceeding 
$ b_{max} $, i.e., having satisfied condition A, we return to the 
discussion of condition B. If B is satisfied, that is, if $ (b_{max} =) 
\: b_{1}' < b_{k-1} $, detrainment of water into layer $ k $ as described 
by (\ref{eq:10}) is the only viable option. If, on the other hand, condition 
B is violated (as illustrated in Fig. 1), it is possible to divide the 
detrained water among two isopycnic layers, $ k $ and $ k-1 $. Detrainment 
can proceed in this case to depth $ L $ without raising the mixed layer 
buoyancy above $ b_{max} $. To implement this, a lower interface depth for 
layer $ k-1, \: z_{k-1}' $, is placed between $ z_{1} $ and $ z_{1}' $ in 
such a way that redistribution of heat in the layer between $ z_{k-1}' $ 
and $ L $ increases the buoyancy in the depth interval $ (z_{k-1}' , 
z_{1}') $ from $ b_{k} $ to $ b_{k-1} $ while decreasing the buoyancy in 
the interval $ (z_{1}' , L) $ from $ b_{max} $ to $ b_{k-1} $. Disregarding 
once again the difference between heat and buoyancy fluxes, the appropriate 
interface depth $ z_{k-1}' $ is found to be
\begin{equation}
z_{k-1}' = z_{1} \frac{b_{1} - b_{k}}{b_{k-1} - b_{k}}
         - L \frac{b_{max} - b_{k-1}}{b_{k-1} - b_{k}}. \label{eq:11}
\end{equation}
In summary, violation of condition B leads to a three-way split of 
mixed-layer water, with water between $ z_{1} $ and $ z_{k-1}' $ being 
detrained into layer $ k $, water between $ z_{k-1}' $ and $ L $ being 
detrained into layer $ k-1 $, and water above depth $ L $ forming the new 
mixed layer with buoyancy $ b_{max} $. The resulting layer configuration 
is indicated by the heavy solid line in Fig. 1.

Momentum exchange directly follows the exchange of mass between the mixed 
layer and the underlying isopycnic layers. The guiding principle is that 
momentum does not vary vertically within a layer, neither before nor after 
vertical redistribution. Total momentum conservation in the column is 
easily achieved under this condition.

A few comments are in order concerning the effect of nonlinearities in the 
equation of state. These nonlinearities are such that transferring a given 
amount of heat from one water parcel in the fluid column to another parcel, 
which may be either identical or lighter than the first one, always
yields a smaller buoyancy decrease in the first parcel and a larger buoyancy 
increase in the second one than a linear buoyancy-temperature relationship 
would suggest. (This is the inverse of the often-noted slight density gain 
associated with mixing two parcels of different buoyancies.)

Thus, a more careful treatment of nonlinearities in the temperature-buoyancy 
relationship during detrainment will have a stabilizing effect on the water 
column if heat is transferred upward, and a destabilizing effect if the 
direction of heat transport is downward. The only place where the latter 
situation occurs in our present detrainment scheme is the (relatively rare) 
three-way detrainment mode, which requires some heat to be tranferred 
downward past the level $ z_{k-1}' $ defined in (\ref{eq:11}). One way to 
avoid potential problems in this case is to evaluate the heating and cooling 
requirements in the slabs above and below $ z_{k-1}' $ independently 
and, realizing that more heat is extracted from above than is used below 
$ z_{k-1}' $, to close the heat budget by adding the positive residual 
to the mixed layer. In summary, no problems should be expected if the 
linear temperature-buoyancy relationship underlying the above detrainment 
formulae is replaced by a more realistic one.

\vspace{2cm}

4. {\bf Atmospheric Forcing Functions}
\nopagebreak

The hybrid model described in the previous sections requires surface 
boundary conditions in the form of time- and space-dependent wind stress 
fields, radiative fluxes and sensible and latent heat fluxes. For the model 
integrations discussed in this paper, a zonally-averaged forcing set was 
constructed for the North Atlantic between $ 10^{0} $ and $ 60^{0} $ N from 
data assembled by Lamb and Bunker (1982) and from the Comprehensive 
Ocean-Atmosphere Data Set (COADS). Because we express the turbulent heat 
fluxes in terms of air-sea temperature differences (Haney, 1971), the 
quantities needed to drive the model include the surface wind speed, the 
surface air temperature, and the near-surface humidity. In order to 
simplify the forcing, we have combined the sensible and latent heat fluxes 
via a procedure described below.

For this analysis, the $ 2^{0} $-square, monthly averaged COADS product was 
used (Woodruff et al., 1987).  We extracted fields of the two horizontal 
wind components, pressure, air and sea temperature, and relative humidity, 
as well as products of these quantities, for the area $ 10^{0} - 60^{0} $ 
N, $ 6^{0} - 74^{0} $ W for the years 1950-1979.  These data were 
time-averaged to extract a 30-year climatological cycle and then zonally 
averaged. In this fashion, monthly, latitude-dependent fields of sea-surface 
temperature $ T_{S} $, zonal wind stress $ \tau_{x} $, wind speed $ U_{10} $, 
and total (latent plus sensible) turbulent heat flux $ H $ were obtained. 
In forming the fluxes from the raw data products, constant transfer 
coefficients for heat, moisture and momentum were used. (The zonally 
averaged meridional wind stress field was also obtained but found to be 
insignificant.) An {\it atmospheric pseudo-temperature} was then calculated 
from
\begin{equation}
T^{*} = T_{S} -\frac  {H} { \rho C_{p} C_{H} U_{10} } \label{eq:12}
\end{equation}
where $ \rho $ and $ c_{p} $ are the air density and specific heat at 
constant pressure, and $ c_{H} $ is the heat flux transfer coefficient. 
$ T^{*} $ is the temperature (considerably lower than the actual air 
temperature) which, when used in the sensible heat flux formula, yields 
a heat flux equal to the combined sensible plus latent heat flux. This 
approach saves the time-consuming computational step of calculating the 
saturated humidity as a function of $ T_{S} $. It implies that, if the 
model calculates a correct zonally averaged sea-surface temperature (SST), 
the surface turbulent fluxes will be identical to the observed fluxes.

While the surface marine observations in COADS provide reasonably good 
estimates of wind, temperature, moisture, and the turbulent fluxes (at 
least for long-term averages), radiative flux data are not included in 
COADS.  Further, while there are a number of atlases that include 
radiative flux data, the accuracy of such data sources, particularly given 
the need for monthly variability, is difficult to assess. Therefore, we 
used the Lamb and Bunker (1982) data to obtain ``radiative fluxes'' 
indirectly, as a residual of the heat budget of the North Atlantic. This 
method has the advantage that unknown modulators of radiative fluxes, 
notably clouds, are implicitly included in the forcing. It further allows 
the basin-averaged energy input to be adjusted to net zero mean. This is 
necessary for the calculations discussed here because the North Atlantic 
exports heat to latitudes higher than $ 60^{0} $ N, and receives heat from 
the tropics south of $ 10^{0} $ N.

In principle, is it possible to project these forcing quantities onto a 
model domain with appropriate interpolation techniques.  However, such an 
interpolation necessarily would be tailored to the time and space grid 
used in a particular integration. Further, the small-scale (i.e., less 
than 3-month or 5-degree) noise in the data is not relevant in the present 
context. Therefore, smoothing and transformation of the data is indicated. 
This was accomplished by developing analytic formulas for the forcing 
using regression techniques. Fields derived from the formulas (plotted 
with the same resolution and smoothing as the original data) are shown 
in Figs. 2-4.

The dynamic forcing by the zonal stress is shown in Fig. 2, with the 
COADS data in Fig. 2a and the patterns from the analytic formula in Fig. 
2b. The formula captures the mid-latitude westerlies, with strongest 
winds in the winter, and less seasonally modulated trades south of about 
$ 30^{0} $ N. Thermodynamically, the model is forced by both the imposed 
radiative flux (Fig. 3) and by the turbulent heat flux. The latter acts, 
in our formulation, as a Newtonian forcing term. (An ``imposed'' flux is 
one that is independent of the model state, whereas a ``Newtonian forcing 
term'' is proportional to the difference between an externally specified 
and a model-generated value. The Newtonian forcing will tend to relax the 
model toward the externally specified value.) The formula for the radiative 
flux yields considerable smoothing of the ``flux'' obtained from the data 
set, but, once again, the main seasonal variability and meridional trend is 
retained.  Likewise, the formulae for $ T^{*} $ and $ U_{10} $, when used 
with the observed SST, give a reasonably good representation of the 
observed turbulent flux (Fig. 4).  In this case, however, there is 
small-scale noise in the result derived from the formulae because the SST 
was used without smoothing.

It is of interest to examine what SST the above forcing would produce in 
the case of a motionless model ocean with zero thermal inertia. This 
equilibrium temperature, $ T_{flux} $, obtained by balancing the turbulent 
and radiative fluxes, is
\begin{equation}
T_{flux} = T^{*} + \frac {R} { \rho C_{p} C_{H} U_{10} } \label{eq:13}
\end{equation}
Fig. 5 compares the COADS SST with $ T_{flux} $. Note that, to the extent 
that the denominators in (12) and (13) can be considered identical, 
$ T_{S} - T_{flux} \propto  H + R $. In physical terms, $ T_{flux} $ is the 
zonally-averaged North Atlantic SST in the absence of entrainment of cold 
water below the mixed layer and without large-scale advection and heat 
storage in the mixed layer. Accordingly, the temperatures in Fig. 5b are 
too warm in the southern part of the domain and too cold in the northern 
part. Note also the agreement in phase between $ T_{flux} $ and the 
radiative heating (Fig. 3b).  Comparing this temperature pattern with the 
model results presented below will show the effects of the modeled currents, 
seasonal heat storage, and vertical entrainment on the temperature field.

\vspace{2cm}

5. {\bf Basic Features of the Modeled Circulation}
\nopagebreak

We report here on the results of a 50-year model integration. In order to 
expedite the spin-up of the basin gyres, the initial density structure 
depicted in Fig. 6 was imposed on a motionless ocean. This zonally 
symmetric structure was excerpted from the Levitus (1982) data set. The 
geostrophic adjustment associated with this density structure, leading to 
a single anticyclonic gyre, occurs on the time scale of $ 1/f $ and is 
completed in a few days. After this initialization, further spin-up occurs 
primarily in response to the wind forcing outlined in the previous section. 
In this section we discuss features of the zonally averaged basin-wide 
circulation, including the annual cycle of mixed layer depth.

A simple, yet effective, way to test whether the model is past the initial 
spin-up stage is to depict a ``census'' of the total mass in each model 
layer as it varies from month to month due to the annual growth and retreat 
of the mixed layer. Such a census for the last 15 years of model time, based 
on a 60-day sampling interval, is shown in Fig. 7. The curves, which 
represent deviations of layer thickness from a long-term mean, indicate that 
the top three layers (i.e., the mixed layer and the two warmest isopycnic 
layers) have settled into a statistically steady state by year 50. The lower 
layers, on the other hand, continue to evolve slowly. We attribute this 
long-term trend to the absense of two processes from the model: 
cross-isopycnic mixing and export of North Atlantic Deep Water into the 
South Atlantic. The gradual transfer of water from layers 4 and 5 to layer 6 
evident in Fig. 7 is the isopycnic model's way of signaling a slow cooling 
trend in the basin. This heat loss is also depicted by the slight downward 
trend in the bottom curve showing annual fluctuations in total heat content. 
The rate at which the model accumulates water in layer 6 (less than 1 m/year 
according to Fig. 7) is inconsequential for the types of processes of 
interest here.

Through a fortuitous choice of isopycnic layer parameters, the wintertime 
mixed layer in the northernmost part of the basin is able to penetrate the 
lowest model layer, thereby enabling the model to produce the model 
equivalent of North Atlantic Deep Water. The rate at which this happens, 
however, has nothing to do with the actual Deep Water formation rate in 
the North Atlantic: neither does the model receive overflow from the
Norwegian Sea nor is it equipped to simulate episodes of subgrid-scale deep 
convection.

The structure of the wind forcing field shown in Fig. 2b implies three 
gyres in our rectangular basin: a subpolar, cyclonic gyre occupying the 
northern part of the basin down to about $ 45^{0} N $ (the latitude of 
maximum eastward wind stress), a subtropical, anticyclonic gyre occupying 
the space between $ 45^{0}N $ and $ 23^{0} N $, and another cyclonic gyre 
south of $ 23^{0} N $. Due to the placement of the southern wall, no 
physical significance should be assigned to the southernmost gyre. Its 
main purpose is to provide a buffer zone between the subtropical gyre and 
the edge of the computational domain.

The strength of the three gyres, indicated in Fig. 8 by the annually 
averaged streamfunction of the barotropic mass flux, is within a few 
percent of the Sverdrup transport implied by the wind forcing. Since the 
model is not eddy-resolving, it can not reproduce the formation of free 
jets and their associated mid-ocean frontal zones. Maximum surface 
velocities in the basin interior are therefore only a few cm/s. The 
western boundary currents on the other hand, which develop strongly 
even in coarse-resolution models because of the need to close the Sverdrup 
flow, reach velocities of up to 30 cm/s in our model.

Because our main objective is the simulation of thermodynamic processes in 
the upper ocean, it is of interest to examine the response of the mixed-layer 
variables in detail. Fig. 9a shows the annual cycle of the zonally 
averaged SST derived from model year 50, to be compared with the observations 
shown in Fig. 5a.  The main differences between the model result and the 
observations is the larger amplitude of the annual cycle in the model. Much 
of this can be traced to the response of the mixed-layer depth shown below. 
The other obvious difference between Figs. 5a and 9a is the sharper 
frontal zone at about $ 40^{0} $N in the observations. We attribute the lack 
of this frontal zone in Fig. 9a primarily to the model's inability to 
generate and maintain free jets equivalent to the North Atlantic Current 
system. Aside from this aspect, the modeled SST is reasonably realistic, 
given the approximations of a rectangular basin with no flow across the 
northern and southern boundaries. In particular, the maximum temperatures 
occur later toward the south in both model and observations, and are quite 
well in phase except in the zone occupied by the truncated tropical gyre. 
Note that the maximum temperatures in the tropics and subtropics occur 
significantly later than the SST maximum inferred from (\ref{eq:13}) 
(Fig. 5b). This shows the effect of the thermal inertia of the upper ocean.

The annual cycle of the zonally averaged model mixed-layer depth is shown 
in Fig. 9b.  As will be noted in the next section, the values are 
reasonable, with the exception that the retreat in the spring occurs rather 
too abruptly. This is due to the smooth nature of the forcing functions 
which do not reflect storm events. In nature, such storms, through wind 
stirring and cooling, produce episodes of deepening during the heating 
season and preclude the monotonic mixed-layer retreat exhibited by the 
model. Despite this overly abrupt retreat, however, the overall annual 
amplitude of mixed-layer depth is quite similar to that observed (e.g., 
Robinson et al., 1979). However, the abrupt retreat implies that the 
mixed-layer's thermal inertia drops too quickly, compared to the real ocean, 
in the early part of the spring heating season; this allows the model SST 
to increase faster than the real SST. 

The annual variability of the zonally averaged model heat content, relative 
to the annual zonal mean, is shown in Fig. 10a. As one would expect, the 
heat intake (during summer) and heat loss (during winter) is most pronounced 
in the northernmost part of the basin. 

Fig. 10b shows the surface turbulent heat flux from the last year of 
the model integration (this is analogous to Fig. 4b except that the 
model SST is used in place of the observations). The main feature that 
is missing from the modeled fluxes is the wintertime maximum at $35^{0} $N. 
This can be traced to the overly large amplitude of the SST annual cycle 
at that latitude (evident, for example, in the behavior of the $ 16^{0} $ 
isotherm).  The lower SST in the model will produce lower heat fluxes, 
and this is particularly evident in winter in the middle of the basin.

In summary, the gross features of the modeled circulation and SST behavior 
are quite reasonable, given the idealizations and limitations of the 
model. The features that depart significantly from observations can be 
traced to specific model compromises, notably limitations in horizontal and 
vertical grid resolution.

\vspace{2cm}

6. {\bf Geographic Variations in the Mixed-Layer Growth Cycle}
\nopagebreak

The seasonal cycle of the oceanic mixed layer depth is controlled by two 
processes: turbulent entrainment due to surface buoyancy and momentum fluxes 
on the one hand, and buoyancy preconditioning due to horizontal advection on 
the other hand. Both processes are included in our model and interact in a 
way that leads, despite the purely zonal nature of the forcing functions, to 
meridional {\em and} zonal variations in mixed layer behavior. 

In this section we will concentrate on the geographic variations in the 
annual water mass exchange between the mixed layer and the interior isopycnic 
layers. We shall focus on conditions in March and in September, the months
that represent the extremes of winter and summer conditions in the layer.
Figs. 11a and 11b illustrate the simulated mixed-layer depth during
these two months. These figures should be compared with the corresponding
charts of observed mixed-layer depth in the North Atlantic, which have 
been published by Robinson {\em et al.} (1979) and by Levitus (1982). 

During March, the bottom of the simulated mixed layer slopes down toward the 
north from about 125m at $ 15^{0} $ to more than 1000m at $ 60^{0} $. This 
agrees reasonably well with the data in the Robinson atlas, which shows 
layer depths of order 50-100m at a latitude corresponding to the southern 
border of our model. At $ 60^{0} $ the Robinson atlas shows extensive 
regions with a layer depth in excess of 900m. The numbers in the Levitus 
atlas indicate a maximum mixed-layer depth of only about 500m in the 
sub-polar regions of the North Atlantic. The discrepancy is hardly 
surprising, considering that the bottom of these deep, cold mixed layers 
is usually not well defined, and considering that neither our model nor 
the atlases take any account of local, transient convective events, which 
can cause occasional extensions of the mixed layer down to the sea bottom 
at more than 4000m depth. Ekman layer convergence in the two grid rows 
next to the northern wall, which causes downwelling in excess of 100m per 
year (see Fig. 12) is undoubtedly a contributing factor.

In September, the mixed layer tends to be less than 50m deep almost 
everywhere, both in the atlases and in our simulations. In the real world, 
the retreat of the mixed layer during summer is accompanied by an almost 
complete suppression of turbulent mixing, particularly during the daylight 
hours of active heating. This transient suppression of turbulence has been 
observed by Chereskin {\em et al.} (1986), Shay and Gregg (1986) and other 
investigators. The deepening of the mixed layer near the western edge of 
the subtropical gyre to 70m is a result of warm water advection by the 
western boundary current. The resulting warm SST anomaly causes a reduction 
in the downward turbulent heat flux which in turn translates into a 
larger Monin-Obukhov length or mixed-layer depth.

The large seasonal excursions in mixed-layer depth have a controlling 
effect upon the modeled motion and temperature fields. The mixed-layer 
velocity and temperature in March and September are shown in Figs. 13a and 
13b. The velocity arrows in these figures represent the resultant of a 
geostrophic and an Ekman component. The relative weight of the two 
components depends on the extent to which the mixed-layer depth exceeds 
$ u^*/f $, the scaling depth of the Ekman layer. It is therefore not 
surprising that the wind stress produces much larger apparent Ekman drift 
velocities in the shallow September mixed layer than in the deep March layer.

The mixed-layer depth fluctuations profoundly affect the horizontal mass 
transport in the mixed layer. As shown in Figs. 14a and 14b, the mixed 
layer in both the subtropical and the subpolar gyre transports about four 
times more water in March than in September, despite the much lower 
flow velocities in Spring. 

The annual entrainment/detrainment cycle can be visualized in a number of 
ways. In Figs. 15 and 16 we show the thickness of layers 2 to 5 in March 
and September. Intercomparison of the four panels in Fig. 16 reveals that 
mixed-layer detrainment in the model tends to form geographically distinct 
pools of water in different isopycnic layers. The controlling variable, 
of course, is the temperature distribution in the mixed layer at the end 
of the winter season. The thickness values for layer 2 indicate that the 
warmest thermocline water in the model is being formed predominantly in the 
western third of the subtropical gyre. While the poor vertical grid 
resolution in the model does not allow us to simulate the formation of 
genuine ``$ 18^{0} $-water", the accumulation of layer-2 water near the 
western boundary undoubtedly occurs for a similar reason.

Poleward advection of warm water by the western boundary current, gradual 
cooling of the surface water on its eastward journey, and a slight southward 
Ekman drift combine to give the springtime isotherms a northwest-southeast 
orientation in the western two-thirds of the basin between $ 35^{0} $ and 
$ 55^{0}N $ (Fig. 14a). Consequently, the bodies of detrained water that 
form in layers 3 and 4 in that region are elongated in northwest-southeast 
direction. Detrainment into layers 2 and 3 takes place in the subtropical 
and tropical gyres (compare Figs. 14a and 16), while sizable amounts of 
water circulating within the subtropical and the subpolar gyre are being 
detrained into layer 4. Detrainment into layer 5 occurs predominantly 
(though not exclusively, as will be seen in the next section) within the 
subpolar gyre. Since the subpolar gyre rotates counterclockwise, the warmest 
water in that gyre is found in the eastern portion of the basin; this is 
where the bulk of layer 4 water is formed. The more frigid remainder of 
the mixed-layer water in the subpolar gyre is being detrained into layer 5. 
The coldest spot in the model ocean is its northwest corner where a small 
amount of water is detrained into layer 6.

By controlling the geographic boundaries of the detrained water masses, 
the mixed-layer isotherms also control the shape of the summertime gyres 
in layers 2, 3, and 4. The similarity between the layer thickness isopleths 
(Fig. 16) and the layer mass-flux streamlines near the outcropping zone is, 
in fact, so close that the streamfunctions will not be reproduced here. It 
is worth noting that the northwest-southeast tilt of the layer mass-flux 
streamlines is far more pronounced in our solutions than in the layer-2 
solutions of Huang and Bryan (1987). Most likely, this difference is due to 
the fact that Huang and Bryan use steady, annually averaged forcing and 
thus do not simulate the sudden subduction of large amounts of water into 
interior isopycnic layers.

The same argument may help us explain the congruence of layer thickness 
isopleths and layer mass-flux streamlines, a feature which seems strangely 
at odds with the subduction flow pattern envisioned by Luyten {\em et al.} 
(1983). The fundamental difference between our results and those of Luyten 
{\em et al.} and of Huang and Bryan is that subduction in our model is 
{\em not} a steady process controlled by Ekman layer convergence, but 
rather an abrupt event taking place only once a year. Our model simulations 
show that the "mode water" lenses produced during this event are not 
continually being replenished from the mixed layer, but exist in relative 
isolation until the mixed layer re-entrains a major part of them during 
the following winter season. Further comments on the difference between 
impulsive and steady subduction, especially its implication on the 
potential vorticity distribution in the seasonal thermocline, can be found 
in Woods (1985). 

A striking view of the annual mixed-layer growth cycle can be gained through 
vertical cross sections. The meridional sections shown in Figs. 17 and 18 
run from the northern to the southern edge through the center of the basin. 
Once again, we show conditions in March and September. As in Fig. 6, we 
display the layer interfaces exactly as represented in the model and use 
some artistic license (i.e., a numerical smoothing process) in drawing the 
isotherms to separate them from the layer interfaces.

We have also added contour lines of a passive ``ventilation tracer" which is 
kept at a constant value in the mixed layer and which follows the mixed-layer 
water into the isopycnic interior during detrainment. Since the mapping 
algorithm is unable to depict sharp gradients and thereby makes it difficult 
to detect the presence of narrow slabs of ventilated water in the isopycnic 
interior, we remove the tracer from the mixed layer in the cross sections 
prior to plotting. As presently implemented, the tracer has an infinite 
lifetime. It was introduced into the model after 10 years of integration 
(and maintained in the mixed layer at 100\% saturation); its distribution 
in the basin at year 50 therefore represents the cumulative effect of 40 
years of thermocline ventilation.

Figs. 17 and 18 are intended to give an impression not only of the annual 
advance and retreat of the mixed layer, but also of the mass exchange 
between the seasonal and the permanent thermocline. The seasonal thermocline 
occupies that part of the cross sectional plane that shows the ventilation 
tracer to be present in September but not in March. The patches of tracer 
remaining in the March section (Fig. 17) reside in the permanent 
thermocline. Some of the details of the tracer distribution cannot be 
understood without knowledge of the horizontal aspects of tracer transport; 
these aspects will be dealt with in the following section.

In an effort to reduce the amount of material displayed in this section, we 
have only shown the mixed layer during its annual extremes. However, we 
feel compelled to add to Figs. 17 and 18 one more cross section illustrating 
the model's behavior during the numerically difficult transition period 
between late winter and summer. As discussed at length in Sec. 3, the 
mismatch between the continuous mixed-layer density field and the discrete 
density structure in the interior leads to periods of relatively slow 
mixed-layer retreat interspersed by detrainment ``bursts", which occur 
during springtime heating whenever the mixed-layer density matches the 
density of one of the interior coordinate layers. In a coarse (vertical) 
resolution experiment such as the one reported here, this discretization 
problem causes ``spikes" in the mixed-layer depth during periods of active 
detrainment. A somewhat extreme example of this phenomenon is given in 
Fig. 19.

Despite the noisy imprint they leave on the mixed layer depth, detrainment 
bursts are not as disruptive numerically as we had originally feared. Since 
they leave the buoyancy structure in the grid column largely intact, they do 
not disturb the velocity field in any major way. Nevertheless, the detrainment 
bursts do not completely ``average out" over the year because the spikes in 
the mixed layer depth immediately get smeared out by the pressure smoothing 
algorithm. Their residual effect on the annually averaged mixed-layer 
entrainment/detrainment rate can be detected between latitudes $ 40^{0} $ 
and $ 53^{0}N $  in Fig. 12. The apparent alignment of the noise elements 
with the isopleths of the late-winter mixed-layer temperature (Fig. 14a) 
and mixed-layer depth (Fig. 11a) is due to the fact that buoyancy profiles 
favoring detrainment bursts are functions of these parameters. 

\vspace{2cm}

7. {\bf Secular Variations of the Tracer Field}
\nopagebreak

The gradual advance over the years of the ventilation tracer in the 
interior isopycnic layers indicates the rate at which the permanent 
thermocline is able to adjust to interannual changes in the thermodynamic 
forcing functions, i.e., to changing climate conditions. The particular 
``climate change" our model is subjected to has to do with the inevitable 
inconsistencies between the observed (and zonally averaged) initial 
conditions and the highly idealized surface fluxes. Until the ventilated 
water, carrying the message about the climate change, has made at least one 
complete journey around the gyre into which it was originally injected, the 
thermocline structure as a whole will continue to evolve -- both in the 
model and in the real ocean. It is this thermocline adjustment process, in 
contrast to the much slower adjustment of the whole ocean basin to a given 
set of climate parameters, which we set out to capture in the present model.

Cox and Bryan (1984) use trajectories to document the rate at which water 
in their model completes one cycle of the thermohaline circulation. In the 
present model, where annual variations in mixed-layer depth are being 
emphasized, the use of such a ``Langrangian" descriptor would be difficult 
because it would require considerable experimentation to locate in the 
seasonal thermocline the narrow slab of detrained water that migrates 
sufficiently far during one summer to remain beneath the mixed layer and 
thereby retain its identity through subsequent winters.

Cox (1985) uses a passively advected tracer (with a finite decay time) to 
describe the thermocline ventilation process in his model. We have, likewise, 
opted for a passively advected tracer as a diagnostic tool. However, by using 
an advection scheme with extremely low horizontal diffusion  (Smolarkiewicz, 
1984), and aided by the fact that diapycnic diffusion in the interior is 
zero, we are able to use a tracer of infinite lifetime. We thereby obtain 
an extremely clear picture of the rate at which ventilated water advances 
in the model ocean. 

The horizontal distribution of the ventilation tracer in the four uppermost 
isopycnic layers 10 and 40 years after its introduction is shown in Figs. 
20 and 21, respectively. In layers 3,4, and 5, the tracer shows clear signs 
of being advected around the subtropical gyre after being injected along its 
northern flank: this is the subduction process treated analytically by 
Luyten {\em et al.} (1983).

Due to the general decrease of current speed with depth, the rate at which 
the tracer circulates around the gyre differs greatly from layer to layer. 
After 10 years, the tongue of ventilated water in layer 3 has advanced south 
and then west near $ 25^{0} $N and is about to be drawn into the western 
boundary current. These results are qualitatively similar to those obtained 
by Cox (1985). No tongue has formed in layer 4 at year 10, but there 
is clear evidence of equatorward migration of ventilated water in the 
eastern third of the gyre. Layer 5 shows no more than a hint of this process. 

After 40 years (Fig. 21), the ventilated water in both layers 3 and 4 has 
made its way around the gyre, encircling pools of unventilated water. In the 
case of layer 3 this pool is rather small and coincides with the position of 
the gyre center in that layer while in layer 4 the pool extends eastward 
from the layer-4 gyre center over almost half the basin width. This 
difference is due to the fact that layer-4 ventilated water is generated 
farther to the east and closer to the northern edge of the subtropical gyre 
than layer-3 water. In layer 5, the tongue of ventilated water in the east 
is barely beginning to turn westward at that time. It appears that the 
water in that layer needs another 75 to 100 years to complete its journey 
around the gyre.

Interpretion of the double maximum in the layer-4 tracer concentration shown 
in Figs. 17 and 18 is now straightforward: the maximum near $ 27^{0} $N is 
due to the westward advection of ventilated water along the southern flank 
of the subtropical gyre while the minimum at $ 36^{0} $N represents the pool 
of unventilated water in the gyre center.

The tracer distribution in layer 2 in Figs. 20 and 21 is relatively 
devoid of features. Water in that layer circulates most quickly and is 
constantly in contact with the mixed layer. The figure does show grid 
regions with less than 100\% tracer concentration, but these are regions 
which during the 40 years of tracer history have never been occupied by 
layer-2 water. (The tracer concentration at each isopycnic grid point 
is initially set to zero, but once ventilated water appears at a grid 
point, the nonzero tracer value remains there even if the water vanishes 
later.) The narrow strip of low tracer concentration next to the western 
boundary in the subtropical gyre in particular attests to the fact that 
layer-2 water is never found on the ``cold" or landward side of the 
boundary current. The low tracer concentration in layer 3 in the same 
location is related to persistent upwelling in the boundary current (see 
Fig. 12).

The layer-6 tracer distribution after 40 years is shown in Fig. 22. Since 
ventilated water in that layer is being generated in the northwestern corner 
of the model basin, the tracer pattern is dominated by signs of cyclonic 
advection around the subpolar gyre. Judging from Fig. 8, the tongue of 
ventilated water seems to extend beyond the boundary of the subpolar gyre. 
However, Fig. 23 indicates that the gyre boundaries in layer 6 (which, 
incidentally, do not show signs of seasonal variability) do not conform 
to those of the barotropic circulation. While the deep water circulation in 
the subtropical and the tropical gyre is aligned with the barotropic gyres, 
and is confined to the narrow recirculation zone near the western boundary 
consistent with the concept of Rhines and Young (1983), the circulation in 
the subpolar gyre not only extends to the eastern boundary, but also occupies 
more latitudinal space than the barotropic gyre. We attribute the dominance 
of the subpolar gyre in the layer-6 circulation to the fact that it is able 
to acquire momentum directly by mass exchange with the wind-forced mixed 
layer. A comparison with Fig. 23 confirms that the tongue of ventilated 
water seen in Fig. 22 has, in fact, not crossed over into the subtropical 
gyre. While this speaks for the numerical precision of the advection scheme 
used, it also indicates that the present model, in failing to resolve 
eddies and thermohaline-driven deep-sea boundary currents, has no means 
of exchanging water across gyre boundaries.

\vspace{2cm}

8. {\bf Summary and Outlook}
\nopagebreak

In this paper we have attempted to show (a) that a Kraus-Turner-type 
mixed-layer model can be successfully grafted upon an ocean circulation model 
framed in isopycnic coordinates; (b) that this combination represents a 
numerically convenient and robust way of subjecting isopycnic models to 
seasonally varying thermohaline forcing; and (c) that a joint model of this 
kind allows us to simulate the process by which fluxes across the air-sea 
interface control the structure and circulation in the thermocline.

In documenting the equations governing mixed-layer entrainment/detrainment 
(in many ways the central theme of this paper) we addressed the problem 
of numerically matching the spatially and temporally variable mixed-layer 
density field to the discrete density structure of the interior ocean. 
We stressed the fact that the detrainment process in a finite-resolution 
isopycnic model, to the best of our knowledge, cannot be implemented in 
thermal-energy conserving fashion without sacrificing either correct mixed 
layer depth or correct mixed layer buoyancy. Inevitably, either choice 
affects, to some extent, the potential vorticity of the detrained water 
masses. 

As with other numerical discretization problems, this one can be suppressed 
by increasing the vertical model resolution. Alternatively, it can be ignored, 
since it more or less remedies itself at the end of each detrainment season. 
Numerous experiments conducted with the model (and its forerunners) over the 
past four years have led us to conclude that low vertical resolution does 
not cause serious numerical stability problems during the detrainment 
process. 

There are a numbers of obvious ways in which the model can be expanded.
Provisions for forcing and advecting the salinity field independently of 
temperature are already in place. (In fact, making sure that the salinity 
retained its initial value of 34.5 mil everywhere in the model domain was 
a convenient way of ascertaining the bulk conservation properties of the 
entrainment/detrainment algorithm.) As the number of prognostic variables 
affecting water density is increased from one to two, an algorithm must be 
added for restoring the density at each interior grid point to its reference 
value. This is done by exchanging small amounts of water among adjacent 
layers, the same process as used by Bleck and Boudra (1981) to restore 
isopycnic conditions in their quasi-isopycnic model. While there appear 
to be no computational problems associated with adding salinity to the 
isopycnic model, for the work reported here we decided against forcing 
the salinity field directly because of the difficulties associated 
with finding the proper forcing functions (particularly the rainfall 
rate).

Our next step, most likely, will be to remove the geographic idealizations 
of the present model and to develop a geographically realistic North Atlantic 
basin model, and to use less idealized forcing. This, in combination with 
the inclusion of salinity as a prognostic variable, will allow the study 
of deep water formation in a more realistic mode, with the potential for 
overspill from the Denmark and Davis Straits and the Norwegian Sea adding 
to the deep water formed by the retreat of the deep wintertime mixed layer. 
This version of the model will be ideally suited to the study of changing 
thermohaline forcing of the oceans associated with climate change (changes 
in freshwater runoff patterns, for example).

As discussed throughout the paper, one limitation of the present formulation 
is the lack of cross-isopycnic mixing. There are no numerical or conceptual 
difficulties surrounding this process which is absolutely necessary (though 
not sufficient in itself) for closing the thermohaline circulation in the 
model. (Horizontal fluxes of abyssal water in and out of the grid domain 
are another major process that will have to addressed in this context.)
In the isopycnic coordinate framework, diapycnic mixing manifests itself 
as mass flux across layer interfaces. Perhaps the quickest way to find the 
mass flux that mirrors the effect of a prescribed vertical exchange 
coefficient $ K $ is to multiply the thermodynamic equation (written here 
in terms of specific volume $ \alpha $) 
\[
\dot \alpha = \left( \frac{g}{\alpha_{0}} \right)^{2} 
\frac{\partial}{\partial p} (K \frac{\partial \alpha}{\partial p})
\]
by $ \partial p / \partial \alpha $. The term appearing on the left-hand 
side of the resulting expression
\[
(\dot \alpha \frac{\partial p}{\partial \alpha}) = 
\left( \frac{g}{\alpha_{0}} \right)^{2} \frac{\partial}{\partial \alpha} 
\left( K ( \frac{\partial p}{\partial \alpha})^{-1} \right)
\]
is the sought-after generalized vertical velocity (in units of pressure 
per time) in $ \alpha $ coordinates. The main challenge, of course, lies 
in the definition of $ K $. 

In addition to the coarse vertical resolution in the present formulation, 
the horizontal resolution used here precluded the formation of eddies 
associated with the Gulf Stream and the North Atlantic Drift. It also limited 
the model's ability to resolve fromtal features observed in the real ocean. 
Future versions of the model will decrease the grid size and be used to 
study these features.

We see the model's primary strength as its ability to portray the basic modes 
of interaction between the oceanic mixed layer and the thermocline even 
under conditions of {\em modest} vertical grid resolution. The experiment 
reported was carried out with only five isopycnic layers beneath the mixed 
layer; yet forcing functions patterned after those observed over the North 
Atlantic produced significant agreement between the numerical solution 
and features of the observed North Atlantic mixed layer if allowance is 
made for the use of a rectangular, flat-bottom ocean basin. The ease at 
which the exchange of mass between the mixed layer, the seasonal 
thermocline, and the permanent thermocline on a seasonal time scale can 
be simulated raises hopes that models of this kind may evolve into viable 
tools for investigating the response of the ocean to climate fluctuations.

\vspace{3cm}

\underline{Acknowledgements:} We are pleased to acknowledge the following 
sources of support: Contract No. N00014-86-C-0363 from the Office of Naval 
Research and Grant No. OCE85-00860 from the National Science Foundation to 
the University of Colorado; and Grant No. OCE86-00593 from the National 
Science Foundation to the University of Miami. Computations were carried out 
at the National Center for Atmospheric Research, operating under NSF 
sponsorship. Among those who contributed to this work through personal 
interaction we wish to thank in particular Dr. John Woods, Hooke Institute, 
Oxford, U.K., and Drs. Claes Rooth and Douglas Boudra, University of Miami.

\end{document}