The effect of different parameterizations and boundary conditions applied to the momentum equation in coarse-resolution thermohaline circulation models.

THIERRY HUCK AND ANDREW J. WEAVER
School of Earth and Ocean Sciences, University of Victoria, P. O. Box 3055, Victoria, BC, V8W 3P6, Canada

ALAIN COLIN DE VERDIERE
Laboratoire de Physique des Océans, Université de Bretagne Occidentale, B. P. 809, 29285 Brest cedex, France

submitted to Journal of Physical Oceanography: Sept. 12, 1996
under the title: The effect of momentum dissipation parameterizations in coarse-resolution thermohaline circulation models.

re-submitted to Journal of Physical Oceanography: June 17, 1997

Corresponding author address: Thierry Huck, Laboratoire de Physique des Océans,
UFR Sciences F3.14, 6 avenue Le Gorgeu, B. P. 809, 29285 Brest cedex, France.
Email: thuck@univ-brest.fr

CONTENT

ABSTRACT

1. INTRODUCTION

2. DESCRIPTION OF THE MODELS
a. The primitive equations
b. The planetary geostrophic equations
c. Choices of momentum dissipation
d. Model geometry, parameters and forcing

3. VALIDITY OF PLANETARY GEOSTROPHIC DYNAMICS

4. GLOBAL DIAGNOSTICS AND THERMOHALINE EFFICIENCY
a. Global, bottom, minimum and zonal averages of temperature
b. Overturning streamfunctions
c. Poleward heat transports
d. Thermohaline efficiency

5. DETAILED COMPARISON OF THE STEADY STATES
a. Horizontal velocity fields and temperature
b. Convection and vertical velocity fields

6. ENERGETICS
a. The potential energy balance
b. Comparison of the energy balance for the different models

7. CONCLUSION

REFERENCES

FIGURE CAPTIONS

TABLE CAPTIONS

ABSTRACT

An ocean circulation model is developed for a Cartesian coordinate flat-bottomed beta-plane, based on the planetary geostrophic (PG) equations, in order to test different parameterizations of the momentum dissipation (Laplacian, biharmonic, Rayleigh and none) and associated boundary conditions. It is used at coarse-resolution for a mid-latitude basin, with restoring boundary conditions for the surface density and with no wind stress. The surface temperature fields and poleward heat transports are quite similar for the equilibrium states obtained using different momentum dissipation parameterizations. However, a comparison of the velocity fields and bottom water properties shows large discrepancies. The traditional Laplacian friction with a no-slip boundary condition produces a satisfying interior circulation, in good agreement with geostrophy and Sverdrup balance, but generates excessively large vertical transports along the lateral boundaries: especially upwelling in the western boundary current (the Veronis effect) and downwelling in the north-east corner. The meridional overturning is thus enhanced, but drives to depth surface waters that are not as cold as the ones in the deep convection regions. Rayleigh friction with a no-normal-flow boundary condition induces a more efficient thermohaline circulation with better agreement between convection regions and areas of downward velocities, colder deep water, much weaker meridional overturning and Veronis effect, but higher poleward heat transport caused by the sharper thermocline. Different closures are used to solve for the tangential velocities along the boundaries but all strongly reduce the (diapycnal) vertical transports in these regions. Although the linear friction lacks physical justifications and is not as satisfying as the Laplacian closure in terms of interior geostrophic and Sverdrup balance, it is an interesting alternative to implement it along with the no-normal-flow boundary condition, since free-slip and no-slip boundary conditions lead to very similar circulations (regardless of the momentum dissipation scheme). Because of the first order geostrophic balance in the ocean interior, the choice of dynamical boundary conditions is more important than the parameterization of the Reynolds stresses for shaping up the thermohaline overturning and the deep water properties in coarse-resolution ocean circulation models: A consistent closure remains to be found to account for lateral boundary layers (with regards to the observations) and is likely to significantly improve models results.

1. Introduction

A fundamental role of the ocean in the climate system is its ability to store and transport heat poleward. In the North Atlantic, the main contribution to the poleward heat transfer is related to the thermohaline overturning, driven by the formation of deep water in the northern seas. Hence, accurately modeling the thermohaline circulation is a crucial task for climate studies. Most current models are run at coarse resolution, the grid-spacing being almost an order of magnitude larger than the Rossby radius of deformation. These models rely on the parameterization of sub-grid-scale processes (e.g. the effect of mesoscale eddies on tracer and momentum dissipation mixing or convection). Although the large scale thermohaline overturning is captured in these coarse-resolution simulations, many aspects of the results need to be improved: thermocline is often too deep; the deep temperature and salinity fields too warm and fresh, respectively and the poleward heat transport in the North Atlantic is typically too weak (Bšning et al. 1995). Veronis (1975) showed that a spurious effect of cross-isopycnal diffusion provides a short-cut to the thermohaline loop by inducing large upwelling of the deep western boundary current into the Gulf Stream. These cold waters that cannot make their way southward significantly reduce the poleward heat transport, hence the whole efficiency of the thermohaline loop is in question. Bšning et al. (1995), however, suggest that this phenomenon will not be avoided unless the Rossby radius of deformation is resolved. Since computing power is not yet available to permit the integration of eddy-resolving models on climatic time scales, another way to reduce this upwelling is searched. A recent parameterization of mesoscale eddy-induced transport of tracers (Gent and McWilliams 1990) seems to help the reduction of the Veronis effect, although this may be limited by numerical constraints on standard advection schemes or computational cost (Weaver and Eby 1997). This parameterization (associated with isopycnal mixing) leads to satisfying improvements concerning the regions of deep-water formation and the water mass characteristics (Danabasoglu et al. 1994; Danabasoglu and McWilliams 1995; Hirst and McDougall 1996). It induces colder deep water (Robitaille and Weaver 1995) due to stronger Antarctic bottom water formation and more pronounced North Atlantic deep water penetration in the Equatorial region, and consequently a less diffuse thermocline. The eddy-induced bolus velocities also reduce the Veronis effect. Veronis (1975) explained the large upwelling at the western boundary in terms of a thermodynamic balance in the tracer equation: large horizontal diffusion across the steeply sloping isopycnals can be compensated only by strong vertical advection, hence breaking the rotational constraint on horizontal divergence. It might alternatively be considered as a numerical artifact of the dynamical boundary condition, since the no-slip implementation within the continuity equation is likely to generate large vertical transports that increase the slope of the isopycnals along the boundaries. This point of view is developed further in this paper (see also Huck et al. 1997). During the last decade, most improvements in global ocean circulation models have focused on the mixing of tracers. On the contrary, the momentum part of these models is scarcely discussed: Greatbatch and Lamb (1990) proposed an interesting parameterization of vertical viscosity such that potential vorticity is mixed along isopycnals. The traditional Laplacian operator copied on molecular viscosity is sometimes replaced by a biharmonic operator for eddy-resolving experiments to increase the scale selectivity of the dissipation (Semtner and Chervin 1988); but at coarse resolution, few alternatives have emerged to parameterize the Reynolds stresses in the momentum equations for representing barotropic and especially baroclinic instabilities (Gent and McWilliams 1996). The departure from geostrophy is also necessary to avoid mass fluxes across the ocean boundaries: the choice between no-slip, free-slip or no-normal-flow boundary conditions are not well justified. In this study, we try to evaluate the consequences of changing the parameterization of momentum dissipation as well as the lateral boundary conditions for coarse resolution thermohaline circulation models. Zhang et al. (1992) demonstrated that a purely geostrophic model with linear momentum dissipation for the barotropic flow can reproduce the results of the Modular Ocean Model (Geophysical Fluid Dynamic Laboratory) for the evolution of tracer fields. In order to clarify this further, a precise comparison of different types of viscosity parameterizations is performed in the same mainframe model, based on the planetary geostrophic equations. To simplify the interpretation, a flat-bottomed Cartesian mid-latitude beta-plane ocean is forced by restoring boundary conditions on the surface density. Wind stress forcing is not included. The steady state solutions are analyzed in terms of large-scale diagnostics, temperature and velocity fields and the global energy balance. The thermodynamics and dynamics explaining the Veronis effect and its sensitivity to horizontal diffusivity, viscosity and resolution is the subject of a later paper (Huck et al. 1997).
The outline of the rest of this paper is as follows: in the next section, we discuss the primitive equation MOM (which serves as our standard reference model) and contrast the assumptions which form the basis of our planetary geostrophic model. We then outline the parameterizations of momentum dissipation utilized for our model intercomparison experiments. Section 3 confirms the validity of planetary geostrophic dynamics while section 4 presents a discussion of the effects of different closure schemes on the oceanic 'climate' diagnostics and the efficiency of the thermohaline circulation under restoring boundary conditions. Section 5 contains a detailed analysis of the temperature and velocity fields. In section 6 we discuss the models' energetics concentrating on the convective-dissipative balance. Finally, section 7 presents a summary and discussion.

2. Description of the models

In this section, we present the simplifications made to the traditional primitive equations (PE) at the origin of most ocean models. For mid-latitude coarse-resolution experiments, we derive the planetary geostrophic (PG) equations, upon which we build a numerical model with various choices for momentum dissipation and associated boundary conditions. We then discuss the choice of the parameters, forcing and methodology for the comparison of the circulations resulting from the different momentum dissipation parameterizations. Table 1 summarizes the different configurations compared in the next section (the letter L' is used for Laplacian viscosity while 'R' refers to linear Rayleigh friction).
The present model intercomparison is carried out in Cartesian coordinates to simplify the set of equations and gain in tractability and understanding. Although the use of these coordinates is not appropriate at the scale of an oceanic basin circulation, we expect only minor deviations in the behavior of the model results when spherical coordinates are employed in our test basin restricted to a mid-latitude beta-plane. As our primary concern is the thermohaline circulation's response to various momentum dissipation closure schemes, the non-linear interaction with the wind-driven circulation is ignored and only thermal forcing is applied.

a. The primitive equations
The formulation of the equations of motion for a Boussinesq, hydrostatic fluid in Cartesian coordinates, using 'large-scale' averaged variables, is: where x (u) refers to the west-east axis (velocity), y (v) to the south-north axis (velocity) and z (w) to the vertical axis (velocity) with positive z, w oriented upward; T is the temperature, r is the density, r0 is the reference sea-water density and a is a constant thermal expansion coefficient (2e-4 K-1 ); f is the Coriolis parameter equal to f0+by; H is the Heaviside step function; zm is the mixed layer depth (assumed uniformly equal to the first level depth); Q is the surface heat flux divided by r0 CP (CP is the specific heat capacity of sea water), parameterized as a restoring boundary condition to a surface climatology following Haney (1971): Q=Q2(T*-TZm/2), where Q2 is a constant related to the restoring time scale, T* is the equivalent restoring atmospheric temperature and TZm/2 is the temperature of the first ocean level. Since we focus on the parameterization of Reynolds stresses in the momentum equations, the mixing of tracers follows the simplest choice, with spatially uniform horizontal and vertical eddy-diffusivities (KH and KV respectively). In order to have a widely used reference and to assess some of the simplifications we make in our planetary geostrophic (PG) model, we compare our results with those obtained with the MOM1 (Pacanowski et al. 1991) for the same Cartesian beta-plane, flat-bottomed, rectangular basin. As in the PG experiments, the wind-stress is not considered, we linearize the equation of state as a function of temperature only and use restoring boundary conditions for the surface temperature. b. The planetary geostrophic equations A scale analysis of the primitive equations for typical coarse-resolution climatic studies shows that time derivatives and non-linear terms can be neglected in the momentum equations; the proper justification for this simplification (Phillips 1963; Hasselmann 1982) requires time scale large compared with a day and space scale large compared with the internal Rossby radius of deformation (~50 km). We can easily assess this simplification with the use of the Rossby number, which measures the ratio of the nonlinear terms over the Coriolis force: Ro=U/(fL)=O(10-2) for typical western boundary velocities of 0.1 m s-1, Coriolis parameter in mid-latitude and grid size of 100 km. Since Robinson and Stommel (1959) and Welander (1959), numerous ocean models have used this simplified set of equations known as planetary geostrophic (PG hereafter) or thermocline equations (see Veronis 1973 for a review):
The tracer conservation equations and equation of state remain the same as in the PE set. Dissipation is usually added in the horizontal momentum equations to represent the sub-grid scale transfers of energy and to allow the velocity field to match the boundary conditions. The relationship between momentum dissipation, dynamical boundary conditions, tracer diffusion and boundary conditions is non-trivial as reported by Samelson and Vallis (1997a), and one needs to be careful to end up with a well-posed set of momentum and tracer equations. The main characteristic of the PG equations is the absence of time derivatives in the momentum equations, which then become diagnostic. Hence, the numerical procedure is straightforward, especially when the barotropic circulation is zero: a baroclinic pressure field is computed from the density field via the hydrostatic equation, the horizontal baroclinic velocity field is then determined from the momentum equations and the vertical velocities are derived from the continuity equation. Only the full prognostic equations for the tracers remain to be integrated in time and the time step is limited by advection only (varying between a day and a week depending on the model dynamics). Zhang et al. (1992) used the PG equations with no explicit friction in the interior (except for the barotropic wind-driven circulation) and no-slip boundary conditions. In comparison to a Bryan-Cox PE model at coarse resolution (under restoring boundary conditions and wind forcing), the simple PG model produced comparable currents and thermocline depth, even though higher vertical diffusivity and lower horizontal resolution were employed in the PE model: they interpreted this as a consequence of the large eddy viscosity in the PE model. Although no friction was added to satisfy the energy transfers from potential to kinetic, the numerical results were reasonable. In the following comparison, we try to remain consistent by using the same resolution, parameters and forcing for each model.
Colin de Verdière (1988, 1989) used this set of equations with the addition of horizontal Laplacian viscosity successfully to check the Sverdrup balance and other traditional theories about the thermally driven circulation; the wind-stress was added later to study the non-linear interaction with the purely thermohaline component of the circulation. This model, coded on an A-grid and solving the horizontal velocity field in spectral space, is used here to validate the B-grid code (consistent with the MOM code) and to evaluate the effects of the numerical grid on the results.
Winton (1993) used linear friction as well as Laplacian friction in similar mid-latitude simulations to examine long time-scale thermohaline oscillations. When horizontal linear friction was used with no-normal-flow boundary conditions, the tangential velocities along the boundaries were determined using the frictional vorticity equation applied to the center of boundary grid-boxes (Winton 1993). This original closure, where the boundary condition applies to the vorticity instead of the momentum equation, is implemented here to provide an alternative to the traditional no-slip and free-slip boundary conditions (both break the rotational constraint on horizontal divergence). Winton and Sarachik (1993) did not analyze in details the discrepancies induced by the momentum dissipation schemes but instead concentrated on climatic time scales; we pursue this work through a detailed comparison of the effects of momentum dissipation and lateral boundary conditions. Salmon (1986, 1990) showed that the use of a linear friction in three dimensions results in a well-posed set of equations if the hydrostatic approximation is relaxed. With a linear equation for the tracers as well, he found analytically the important "physical features" of traditional models (coastal boundary layers), although the lateral non-hydrostatic boundary layers were too narrow to be resolved explicitly. He showed some thermocline simulations in his latter paper, but never applied realistic forcing over a whole ocean basin. Here, we implement this non-hydrostatic set of equations and observe a tendency towards convergence to Winton's closure when the friction coefficient in the vertical momentum equation approaches to zero (exp. 15, 16 and 17 in Table 2).
In order to compare these momentum dissipation parameterization, we consider a modular program iterating the tracer evolution by advection, diffusion and convection, on a B-grid, using the traditional energy-conserving schemes in relation with the hydrostatic pressure formulation (Bryan 1969; Semtner 1986). The vertical velocity field is computed from the continuity equation at the interface of the model levels - pressure and temperature being defined at the mid-levels in the centers of the model grid boxes and horizontal velocities at the mid-levels corners. Depending on the type of momentum viscosity and lateral boundary conditions, we compute the baroclinic horizontal velocity field from the hydrostatic baroclinic pressure field. No wind stress, topography and bottom friction makes our barotropic streamfunction vanish, such that the baroclinic fields define the whole solution.

c. Choices of momentum dissipation
We describe here the method used to solve for the different choices of momentum dissipation, starting with the models most similar to the MOM and proceeding to the most dissimilar. PGL (Planetary Geostrophic model with Laplacian friction) is based on the same equations as Colin de Verdire (1988, 1989), but written on a B-grid and employing the iterative method for solving the horizontal velocities proposed by Winton (1993). From the equation for horizontal momentum, we derive a single complex equation for (u+i v).
We solve this 2-dimensional elliptic equation at each model level by Successive Over-Relaxation for the interior velocities, the boundary speeds being set to zero in application of the no-slip boundary condition. This is not as computationally efficient as Colin de Verdière's direct solution in spectral space, but much simpler to implement (it can also easily handle irregular geometry). Furthermore, it permits lower viscosity coefficients without generating numerical noise. Laplacian friction is widely used in numerical modeling since its scale selectivity damps out small-scale perturbations. Similar to molecular viscosity in its stress-formulation of surface forces, its mixing scheme conserves energy in the interior and dissipates it along the boundaries by lateral stress.
PGLA (Planetary Geostrophic model with Laplacian friction on A-grid) is the model used by Colin de Verdière (1988, 1989). All the variables (u, v, w, T) are computed on the horizontal A-grid and the horizontal velocities are computed by direct inversion of the tridiagonal matrix resulting from the fast Fourier transform of the spatial fields. No-slip boundary conditions are used, such that only the horizontal grid for horizontal velocities differs from the previous model PGL. The scaling of the Munk boundary layer is straightforward from the vorticity equation: . Such a scaling leads to a dissipation length scale of (AH/beta)^1/3. We follow Colin de Verdière's criterion to determine the viscosity coefficient as a function of the resolution Dx: AH=[1.6-2]*Dx^3=1.5e5 m2 s-1 for Dx=160 km and b at 40°N.
PGLslip (Planetary Geostrophic model with Laplacian friction and free-slip boundary condition) uses the same momentum dissipation as PGL, but with free-slip boundary conditions: u.n=0 and d(u.k)/dn=0 on the lateral walls, where n and k are vectors respectively normal and tangential to the boundary. Stewart (1989) suggests that both boundary conditions are equally justified at the resolution we use but points out that the no-slip choice results in no net vorticity and meridional vorticity transport in the western boundary current. Dynamics are solved by the same method as PGL, with the velocities tangential to the boundary equal to the value at the closest inner point for each SOR iteration.
PGBslip (Planetary Geostrophic model with Biharmonic friction) is inspired by high resolution "eddy-resolving" models, as it uses a biharmonic operator for the horizontal friction, associated with free-slip boundary conditions. It is even more scale-selective than a Laplacian formulation, but is known to generate unfortunate local extrema of the velocities. The equations for the horizontal momentum become: .
From the relative vorticity (z=vx-uy) equation: bv-fwz = -A4 (zxxxx+zyyyy), the dissipative boundary layer can be scaled as (A4/b)1/5. Its resolution by our grid-spacing defines an otherwise unphysical friction coefficient: A4 > beta*(2Dx)^5 = 2e-11 m-1 s-1 * [2*160 km]^5 = 1e17 m4 s-1. This value produced a too narrow boundary current and was increased to 1e19 m4 s-1 for the comparison. PG0 (Planetary Geostrophic model with no friction) implements the purely geostrophic circulation with no slip boundary conditions on lateral walls. Therefore, in the interior of the domain, the horizontal velocities u and v are computed respectively as -py/(fr0) and px/(fr0); they are both set to zero on the lateral boundaries. From the continuity equation, we can expect to satisfy the Sverdrup balance in the interior, but also to create large vertical velocities along the boundaries, of order UDz/Dx=O(1e-5 m s-1). These vertical velocities usually reduce the allowable time step to one day in our model configuration. Zhang et al. (1992) point out that this parameterization depends on the grid-spacing and no convergence of the results should be expected when the resolution is refined. This is an important "philosophical" discrepancy with the previous models, whose momentum mixing formulation makes convergence possible. Formally, this formulation is not "well-posed" mathematically but the introduction of numerical boundary layers due to the imposition of the no-slip condition at the coast leads to realistic solutions. This model is also used in conjunction with a free-slip boundary condition (PG0slip) and a no-normal-flow boundary condition as described further with a non-zero linear friction just next to the boundaries (PG0W). PGR0 (Planetary Geostrophic model with Rayleigh friction and "no-slip" boundary conditions) implements linear horizontal friction with zero velocities imposed at the lateral boundaries. Then, the velocity field is simply a linear combination of the horizontal derivatives of the baroclinic hydrostatic pressure: . As in the Zhang et al. (1992) model, the velocities on the boundaries are simply set to zero. From the vorticity equation: , we get a scaling of the width of the Stommel western boundary layer as eH/b. If we want to properly resolve this boundary layer, our horizontal grid spacing needs to be greater than this typical scale. Conversely, for a given resolution of 160 km, eH needs to be greater than bDx (~3´10-6 s-1). Assuming this limiting value, the non-dimensional ratio eH/2W (where W is the rotation rate of the earth) equals 0.03, hence the geostrophic balance holds within a few percent in the whole basin. Like PG0, this model assumes a resolution-dependent closure that lacks the physics to match the boundary conditions. It is also implemented with a free-slip boundary condition (PGRslip).
PGRW (Planetary Geostrophic model with Rayleigh friction and Winton's vorticity closure) uses the same friction as PGR0 with a no-normal-flow boundary condition. The new unknowns are then the tangential velocities along the lateral walls, which we determine with the method introduced by Winton (1993). For each level, we write the frictional vorticity equations at the center of each grid-box along the boundaries, using the unknown tangential velocities. This provides as many equations as unknown velocities for each level and just requires solving a wrap-around bi-diagonal matrix. The vorticity equation is the same as in PGR0. Writing these equations for the center of the tracer boxes (since the velocities are defined at the corners) requires an average of velocities and derivatives carried at the corners and sides of each model box, respectively. An important aspect of this formulation is the necessity of friction for the well-posedness of the problem (the determinant of the system for boundary velocities cancels if eH=0). Moreover, this formulation is not successful for solving a cross-equatorial basin circulation. This procedure generates large horizontal velocities on the boundaries, but reduces the vertical velocities along the vertical walls compared to the no-slip boundary conditions. The scaling of this frictional vorticity equation gives the diagnostic vertical velocities to compensate for the horizontal divergence: compared to 1e-5 m s-1 with the no-slip boundary conditions. Thus, a flow hitting a lateral wall recirculates horizontally in the boundary current, instead of vertically in an intense downwelling or upwelling when no tangential velocities are permitted. Since the tangential velocities increase when the grid-spacing decreases, realistic results require that the resolution remains coarse.
PGR4 (Planetary Geostrophic model with Rayleigh friction and biharmonic tracer diffusion) uses linear friction as well with a no-normal-flow boundary condition, but the tangential velocities along the boundaries are determined through the addition of a biharmonic operator for the tracer diffusion. For the details of this astute formulation, the reader is referred to Samelson and Vallis (1997ab): The biharmonic diffusion coefficient (K4=4.4e14 m4/s) is chosen according to these references to reduce the diapycnal fluxes by horizontal diffusion in the western boundary current. Finally, PGS (Planetary Geostrophic model with Salmon's 3-d linear friction) is based on the full linear equations used by Salmon (1986, 1990):
Salmon (1990) and Samelson and Vallis (1997ab) resolved this system in numerical models with eH=eV. Salmon (1986) worked out a boundary layer analysis where he pointed out the classical Stommel western boundary, of scale (eH/b)=O(200 km) and a non-hydrostatic lateral zone, of scale (eVH/f)=O(100 m). Willing to resolve the latter without refining our coarse resolution, we decided to use a friction coefficient greater in the vertical than in the horizontal: we managed to increase the latter length scale such that it would be partially resolved by our resolution, without perturbing significantly the hydrostatic approximation in the interior of the basin. The reason for the artificial vertical friction is the well-posedness of the system: maintaining hydrostasy as well as horizontal frictional-geostrophy would entirely determine the horizontal baroclinic velocity field, which could never match the no-normal-flow boundary condition from the bottom to the surface of the lateral walls. An alternative solution is used in the previous model (Samelson and Vallis 1997ab) to avoid the resolution of a 3-dimensional problem: retaining the hydrostatic approximation but adding a biharmonic diffusion of temperature to allow no heat flux across the boundaries.
Nevertheless, relaxing the hydrostatic approximation is very costly numerically, as one shifts from a 2-d (per level) to a fully 3-d problem. The two ways of defining the elliptic equation for pressure are to write the velocities as functions of the first order derivatives of pressure and density from the momentum equations, and use these expressions either in the continuity equation or in the vorticity equation. This leads to the following elliptic PDE: The no-normal-flow boundary conditions are given in terms of the normal and right-handed pressure derivatives (respectively noted by the subscripts n and s) from solving the horizontal momentum equations for the horizontal velocities (similar to PGR0): .
These are oblique boundary conditions where the tangential derivative is weighted by a much larger coefficient than the normal one, a characteristic that makes solution even more complicated as it generates large off-diagonal coefficients in the resolution matrix. The idealized geometry allows a simple decomposition of the solution in vertical modes for pressure. The equations derived for each vertical mode and the details of the solution technique can be found in Huck (1997).

d. Model geometry, parameters and forcing
The Cartesian beta-plane, centered at 40°N, extends over 4480 km latitudinaly from 20°N to 60°N and over 5120 km meridionaly (roughly 60¡ in longitude at 40°N). The horizontal resolution is 160 km and the depth of the basin (H) is 4500 m, discretized by 15 levels (respectively 50*3, 100,150, 200, 250, 300, 350, 400, 450, 500, 550*3 m deep). Only PGS differs with 6 times more vertical levels of uniform 50 m depth. The tracer mixing uses traditional constant coefficients for the resolution, with horizontal eddy diffusivity KH=700 m2 s-1 and vertical eddy diffusivity KV=1e-4 m2 s-1. Convection is handled explicitly at each time step by the Rahmstorf (1993) scheme. The momentum dissipation coefficients are derived from the scaling of the boundary layers for each parameterization: AH=1.5e5 m2 s-1 with Laplacian friction and eH=3e-6 s-1 with linear friction, such that eH/2W=0.03. For PGS, different values of the vertical friction coefficient eV are tested (88.5 and 8.85 s-1, equivalent to 0.3 and 0.03 in the scaling of the hydrostatic equation: eVscaled=eVH2/(2Wa2), where a is the earth's radius).
The atmospheric forcing consists of restoring the sea surface temperature to zonally-uniform reference temperatures, linearly varying from 25°C at 20°N to 2°C at 60°N. The restoring time scale is based on Haney (1971) estimated fluxes of 35 W m-2 K-1 (equivalent to a 66 days restoring time constant for the 50 m deep surface layer). Additional tests with fixed flux (Neumann) and fixed surface temperature (Dirichlet) boundary conditions indicate our conclusions on the effect of momentum dissipation are only slightly affected by this arbitrary choice. We chose to present the results from the experiments with the 66 days restoring time scale because it is the most common surface boundary condition for ocean models: these experiments lead to a steady-state solution with realistic surface fluxes (zonally-uniform heat fluxes usually drive decadal oscillations, while Dirichlet condition induces unrealistically high fluxes in the Gulf Stream region).
The time stepping is a leapfrog scheme for the advection and Euler backward scheme for the diffusion, occasionally mixed by forward Euler time steps (Euler-backward in the MOM) once every 7 to 17 time steps, except in PGLA where a Matsuno time-stepping proved faster, although it requires the evaluation of diffusion and advection terms twice every time-step. The allowable time step depends largely on the momentum dissipation choice and boundary conditions, extending from 1 day with the MOM code or PG0 to 6 days in PGS or PGRW. The limiting factor is often the CFL criteria for the vertical advection (but for the northward advection in PGS/PGRW) in the PG models, that has been related to the speed of frictional waves along the boundaries (Colin de Verdière 1988; Killworth 1985). The baroclinic instability in PG equations described by Colin de Verdire (1986) did not occur in our experiments, because of the strong surface restoring boundary conditions and the diffusion of tracer and momentum. Starting from an initial uniform temperature of 4°C in the basin, we ran the models to a steady state (3000 years), by which time the residual surface heat flux varies between 10-3 and 10-5 W m-2. PGL, PGLA and the MOM codes are virtually identical since they implement the same dissipation scheme. Through the analysis of global diagnostics as well as more detailed patterns, we found a large dependency of the results on the dynamical boundary conditions. With reduced vertical transports along the boundaries, PGRW, PGR4 and PGS are quite similar compared to all the other models. The comparison of the results from these two groups leads the discussion: one represents the traditional PE models at coarse resolution with no slip boundary conditions, the other being an alternative emphasizing horizontal instead of vertical recirculation along the boundaries, via tangential velocities along the lateral walls.

3. Validity of planetary geostrophic dynamics

In order to quantify the effect of the terms we neglect in the simplified diagnostic equations for momentum, we conducted a few experiments with the full primitive equations of the MOM. First, we compare an experiment with typical vertical viscosities used in large scale numerical models (~1e-3 m2 s-1) and one with no vertical viscosity. This leads to mean differences of 6e-4 °C in the steady state temperature fields, no visible effect on the transient evolution of kinetic energy during the spin-up and a change in many diagnostic values smaller than 0.1%, as shown in the first lines of Table 2. This is consistent with a rough scaling of the ratio of vertical over horizontal viscosity, assuming the minimum horizontal viscosity imposed by the numerical resolution of the western boundary current (physical values would be much lower though): . The choice of convection scheme is known to have significant effects on the thermohaline circulation (Marotzke 1991; Rahmstorf 1994): to avoid any discrepancies in our comparison, we use the same scheme in all the models. To determine the order of magnitude of its influence, we run the MOM with the implicit vertical mixing scheme and with the full convection scheme written by Rahmstorf (1993). The mean change in temperature is found to be about 2e-3 °C and other diagnostic quantities changed by less than 1%. Note that this is already more than the effect of vertical viscosity, emphasizing how important the numerical representation of convection is. The Rahmstorf scheme is retained in later experiments, as it totally removes static instabilities at the lowest numerical cost (Pacanowski 1995); with the implicit vertical mixing scheme, the temperature profiles in the convection regions were still increasing by 10-4 °C from the surface to the bottom.
At low Rossby number [Ro=U/(f L)=0.1 m s-1/(10-4 s-1(100 km)=1e-2], we expect the non-linear terms in the momentum equations to be negligible in comparison with the geostrophic balance over the whole mid-latitude basin. This assumption, justifying the PG system, is the most important one to check in comparison with the influence of the momentum dissipation scheme. Once again, the discrepancies agree with the scaling: considering the coding and complexity differences between the MOM code and our much simpler PGL model, they produce steady-state solutions that compare very well (a mean temperature difference of 2e-2 °C and a global change around 1% in all the diagnostics). Our PG ocean is slightly warmer and more energetic, it transports more heat poleward and more water in the western boundary current. However, the transient behavior with the full momentum equation is not exactly similar to the PG behavior, especially during the very energetic 30 years of spin-up.
The most surprising comparison was with the Colin de Verdire model, using the same equations written on an A-grid (horizontal velocities and temperature are both in the center of the grid-boxes). Globally, the difference from the B-grid to the A-grid formulation are more important than the effect of the non-linear terms (2.5e-2 °C in the mean temperature, with some diagnostics closer to the MOM and others to PGL).
While the effects of non-linear terms in our experiments is in the range of numerical differences in the coding, which agrees with the conclusions of Weaver and Sarachik (1991), it is relevant to restrict this statement to mid-latitude coarse-resolution simulations.

4. Global diagnostics and thermohaline efficiency

a. Global, bottom, minimum and zonal averages of temperature
Mean basin, mean bottom and minimum temperatures are highly correlated in all the models (see Table 2): in steady state, the cancellation of the area average of the surface heat flux requires that the mean surface temperature equals the mean restoring temperature, i.e. 13.5 °C for these experiments. The coldest temperatures reached at the surface slowly fills up the whole deep basin through convection, followed by deep advection, driving in this way the bottom and basin means. The coldest temperatures occur along the northern boundaries (where the restoring temperature is the coolest), usually in the western half of the basin, as a balance dominated by north-south diffusion, surface forcing and convection. We observed that the horizontal diffusivity is an important parameter influencing the minimum surface temperature in these models, as seen in Fig. 1: a striking linear relationship exists between the mean basin temperature (or bottom temperature) and the horizontal diffusivity in a log-log plot, especially for PGL. For PGRW, two regimes with different slopes seem to occur below and above 700 m2/s, the latter one being very similar to PGL. As concluded by Bryan (1986), the sensitivity to KH is due to its effect on cross-isopycnal mixing, and as such, an artifact of the horizontal instead of isopycnal formulation. We implemented isopycnal mixing in the MOM to evaluate the implied modifications on the temperature field: The bottom temperature cools to 3.16°C (3.53°C originally for the MOM) but a background horizontal diffusion of 100 m2/S is required for numerical stability. The addition of the Gent and McWilliams (1990) parameterization for mesoscale eddy-induced tracer transport (isopycnal thickness diffusion of 1000 m2/s) allows to eliminate the background diffusivity and the mean bottom temperature cools further to 2.92°C (3.24°C with horizontal mixing instead of isopycnal mixing).
We can now highlight the significant differences achieved by the different models (for the same horizontal diffusion) in terms of bottom temperature. PGRW and PGS produce the lowest values (around 3.4 °C), PG0 the highest (3.75 °C), the Laplacian models being in the middle (~3.54 °C). The mean basin temperature varies consequently by 0.5 ¡C depending on the momentum dissipation and boundary conditions used. Such a variation is more than the effect of the change in KH from 500 to 1000 m2 s-1 within the same model (0.16 ¡C in PGRW and 0.25 °C in PGL). These significant differences justify a posteriori our study: viscosity parameterization and lateral boundary conditions have a strong influence on the tracer field, comparable to the effect of varying the tracer mixing coefficient within its range of uncertainty. Among these two dynamical choices, the change in the lateral boundary condition appears to be the main factor, as the interior dissipation accounts for only 0.05 °C from PGR0 to PGL, as long as it is not purely geostrophic (PG0 differs from these by 0.3 °C). We delay the interpretation of these results until we link them with the dynamics. Additional factors influencing the basin mean temperature (via minimum and bottom temperature) are the horizontal resolution and dissipation coefficients, because of the local balance between horizontal advection-diffusion, surface cooling and convection leading to the coolest surface waters. Increasing the horizontal diffusivity warms the minimum surface temperature (by smoothing the surface temperatures). Accordingly, decreasing momentum dissipation increases kinetic energy, hence the magnitude of the advective terms in the tracer balance compared to the surface fluxes: the minimum temperature increases, although the relaxation cooling increases as well. Increasing resolution (hence the role of advection relative to diffusion) has a similar effect.
Figure 2 shows latitude-depth sections of the zonally averaged temperature in the thermocline. Straub (1996) points out that this field is not very sensitive to the internal dynamics of the model, since the boundary currents concern only a small ratio of the basin area. The main effect we attribute to the dynamics is to set the bottom water temperature and hence the temperature at the base of the thermocline, since the deep waters are very weakly stratified. The linear friction model with no-normal-flow boundary condition induces the sharpest thermocline, but in disagreement with the conclusions of Zhang et al. (1992), the thermocline is more diffuse in the purely geostrophic model than in PGL/MOM (this is even more clearly represented by the zonally averaged Brunt-VŠisŠlŠ frequency). We do not observe the significant difference in the meridional overturning between PG0 and PGL that support their explanation, although it is supported in terms of cross-isopycnal transports (the maxima of the cross-isopycnal overturning streamfunction are respectively 11.2 and 13.1 Sv for PGL and PG0). Since their conclusion holds with no wind forcing, we must consider the influence of the horizontal grid (D-grid in their model) or the details of the coding (S. Zhang personal communication). The upwelling along the western boundary is twice stronger in PG0 than in PGL, even in terms of diapycnal upwelling. Then, the more diffuse thermocline in PG0 appears as a consequence of intensified vertical mixing by upward and downward motions and is correlated with the bottom water temperature. The narrower extent of the convection regions in PGRW certainly plays a role in the efficiency to produce colder deep water, but the vertical velocity fields achieved by the various models also modify the global equilibrium of the thermocline. Since the vertical transports occurring along the boundaries are of the same order of magnitude as the ones in the interior (section 5b), the lateral boundary condition is more important than the momentum dissipation, which modifies slightly the horizontal divergence in the interior.

b. Overturning streamfunctions
The meridional overturning streamfunction (MOSF) is defined as a vertical integration of the zonally-averaged northward velocities: .
However, we should be aware of the purely dynamical aspect of such a streamfunction, which gives an insight of the mass transport, but by no way of the heat transport. The zonal overturning streamfunction (ZOSF) is analogously defined as a vertical integration of the meridionally-averaged eastward velocities: . These purely dynamical diagnostics are certainly among the most sensitive to the type of viscosity in our experiments, as seen in Figs. 3 and 4. The MOSF varies from 7.6 Sv (1 Sv=106 m3 s-1) in PGS to 15.6 Sv for PG0; once again, these models draw attention to the most extreme variations in Table 2. The mass transports show the greatest sensitivity to the type of boundary conditions, not to the momentum dissipation. The no-normal-flow boundary condition with vorticity closure generates a MOSF below 9 Sv, the free-slip one with harmonic or biharmonic viscosity, around 11 Sv, and the no-slip one, above 14 Sv (whatever internal dissipation is chosen). Just the closure of lateral velocities along the boundaries, whose divergence feeds the vertical velocity field, can therefore largely modify the globally averaged diagnostics, by as much as 100%.
We observe the same variability between experiments with the zonal overturning streamfunction, showing generally the same pattern of two loops, one concentrated at the surface and extending along the western boundary, the other deeper along the eastern boundary. Both are associated with upwelling along the meridional borders, especially the maximum of the zonal overturning is almost entirely set by the maximum upwelling along the western boundary. These two loops differ significantly in intensity and spatial distribution according to the boundary conditions. This is shown in Fig. 4 where the eastern deep cell transports from 2 Sv in PGRW and fills the whole width of the basin, to 4 Sv in PGL where it is deeper and restricted to less than 2000 km longitudinally. In agreement with the thermal wind balance, these patterns are coherent with the north-south density gradient, which is negative only in the regions where the ZOSF is negative. The surface loop is even more variable, with 2 Sv in PGS/PGRW (near the surface) and up to 8 Sv in PGL(A)/MOM (extending to the bottom of the western half of the basin), again showing the preponderance of the boundary closure, especially along the western boundary.
We estimated the northward transport achieved by the western boundary current in the different models and the latitude of its maximum (columns 9 and 10 in Table 2). The method involves finding for each latitude the maximum of the northward transport, integrated from the western boundary and the surface, as a function of longitude and depth (except in PG0, it is often similar to the northward transport across the upper 850 m ´ western 640 km area). We found a good correlation of this value with the MOSF in most of the cases: this confirms the traditional thermohaline loop where the WBC brings the water that will sink in the northern regions and feeds the deep WBC. Therefore, the separation of the thermocline circulation into its divergent and rotational parts shows a huge preponderance of the former. In this respect, the lack of wind-forcing and non-linear terms in the momentum equation is particularly important. c. Poleward heat transports At equilibrium, we can rewrite the vertically- and zonally- averaged form of the tracer equation showing the balance between the advective-diffusive heat transport and the surface heat flux. Integrating from the southern boundary (where no flux occurs) to any latitude, we can define the advective and diffusive poleward heat transports, whose sum is the total ocean heat transport (equal to the integrated surface heat fluxes at equilibrium): The maximum of this value measures the efficiency of the ocean to transport heat northward, for a given surface forcing. This value seems more dependent on the restoring time scale than on the internal dynamics of the ocean. Choosing the horizontal diffusivity (representing the non-resolved eddy-induced heat transport) then determines the diffusive part of the heat transport and the advective part is that which remains. The vertical diffusivity has a major influence on the overturning circulation and the poleward heat transport, by setting the thermocline depth and providing the potential energy to the whole system (Bryan 1987; Colin de Verdire 1988; Bryan 1991). To avoid any effect of this critical tracer mixing coefficient, we kept it uniform and constant throughout the experiments. Finally, the range of variation of the maximum advective heat transport remains around 15%, from 0.21 PW (1 petawatt=1015 watt) for PG0 to 0.24 for PGRW/PGS (Table 2). The atmospheric feedback of the ocean is then weakly sensitive to the momentum dissipation parameterization, at least compared to the meridional overturning.

d. Thermohaline efficiency
A fundamental role of the ocean in regulating climate is its ability to store and transport heat poleward. In the mid-latitude Atlantic, the overturning transport in the meridional-vertical plane dominates (Bryan 1962). The amount of heat released to the atmosphere by this circulation is approximately given by the product of the overturning rate and the temperature change required to convert poleward flowing water to deep water:
We therefore wish to consider the efficiency of the PHT as compared to MOSF and bottom water temperatures. Zonally-averaging our ocean, we find a thermocline with northward velocities generally in the first 800 m and deep water with southward velocities below, as observed in the western boundary. The PHT of such a simple 2-layer model is proportional to the transport and the temperature difference between the upper and the lower branch. For an estimation of DT, we can use the difference between the mean thermocline temperature (around 9°C since the mean surface temperature in steady state equals the mean restoring surface temperature, i.e. 13.5°C) and the mean bottom temperature (around 4°C). This gives a rough scaling of the expected PHT: 4e6 J m-3 K-1 ( 5°C ( 10 Sv = 0.2 PW (1015 W). To assess this relationship, we examine the maximum advective PHT as a function of MOSF and DT (varying as the opposite of the mean bottom temperature) for the different models. Although the correlation with MOSF (Fig. 5a) is better than with the bottom temperature (Fig. 5b), both are negative. From the last formula, we were expecting the opposite correlation between PHT and MOSF. But considering the strong correlation between MOSF and mean bottom temperature for all the models (Fig. 5c), the shallow to deep water temperature contrast controls the PHT in spite of the change in the overturning rate. In analogy with electricity transport in wires, the higher the voltage, the lower the loss of energy. And it is the overall efficiency of the thermohaline loop that enables the ocean to transport heat poleward However, these correlations may be a consequence of the restoring boundary conditions, since the few experiments we did with zonally uniform surface heat flux seem to produce a negative correlation between MOSF and mean bottom temperature (r=-0.81 for 10 different models, giving a low 95% confidence level of 0.69). In these heat-flux boundary condition experiments, the relationship between the maximum advective PHT and MOSF or DT was consistent with the restoring experiments, although the correlation coefficients were lower than the confidence level. Since the mean surface temperature is no longer fixed with this forcing, it is difficult to interpret these results. Under restoring boundary conditions, there is a strong constraint on the surface temperature field, but the advective PHT is more free: when the deep water temperature is warmer, then the overturning rate gets stronger; this process can change the PHT. With constant heat-flux boundary conditions, the constraint is on the total poleward heat transport (on the advective part as well since it is the major contributor). Depending on the dynamic part of the model, the surface temperature can change along with the bottom water temperature to produce the right heat transport. Hence, we observe a lower variability of the mean bottom water temperature.
These results suggest that the thermocline to deep water temperature contrast is more important than the overturning intensity in controlling the poleward heat transport. Both should be considered when assessing a global model's efficiency to transport heat poleward. Finally, we look at the correlation between the western upwelling (as a measure of the Veronis effect) and the deep water temperature (Fig. 5d). The extremely high correlation coefficient clearly shows that the short-cut of the thermohaline loop (via upwelling in the Gulf Stream around 50 °N) does primarily influence the deep water temperature and consequently the whole thermohaline efficiency (Bšning et al. 1995).

5. Detailed comparison of the steady states

a. Horizontal velocity fields and temperature
Figure 6 shows the steady-state horizontal circulation and temperature averaged over the thermocline depth (determined as the mean depth of horizontal current inversion, around 850 m) for PGRW, PGL and PG0. Although the main gyre pattern is similar, there are some noticeable differences. The western boundary current has a much broader and more diffuse signature in PGRW than in PGL or PG0. This is correlated with the typical recirculation east of the WBC. Associated with this southward mass transport, there is also a significant southward heat transport. The overshoot of the Gulf Stream is characteristic of the Laplacian dissipation (Munk 1950), and also of biharmonic friction. The linear friction solutions lack this feature, and produce a weaker northward transport of mass, consistent with the differences in the barotropic wind-driven circulation between Stommel (1948) and Munk (1950). The tropical thermocline reaches higher temperature with PGL and PG0 although the horizontal circulation in this area is stronger than in PGRW. In the north-east corner, the temperature field and the velocities differ significantly and illustrate the main effect of the boundary velocities associated with the vorticity closure. Over the eastern half of the northern basin, PGL and PG0 isotherms are perfectly zonal so that the eastward current hits the wall perpendicularly and descends to depth in an intense downwelling. On the contrary, PGRW shows northward curving isotherms and currents in the north-east corner: the impinging eastward current recirculates in the first 250 m in northward tangential velocities along the eastern wall, but this recirculation generates southward boundary currents below and the depth-average over the thermocline almost entirely cancels. The westward velocities along the northern boundary in PGRW, PGR4 and PGS are rather unexpected, but they contribute to the recirculation initialized in the N-E corner. Winton and Sarachik (1993) already commented on these, and found them consistent with the observations of McCartney and Talley (1982) under a more realistic forcing. However, realistic topography and coast line would certainly profoundly change this feature by adding the varying depth in the vorticity constraint. These boundary velocities are responsible for the westward shift of the coldest regions from PGL to PGRW as described below. At depth, PGL has a much larger range of temperatures than PGRW (0.036 °C standard deviation compared to 0.008 for the 4225 m deep level), the waters driven by the downwelling in the north-east corner being 0.2 °C warmer than the ones formed by convection along the northern boundary, close to the north-west corner. However, the horizontal circulations at depth are fairly similar (except in the N-E corner, and stronger in PGL), since they are primarily driven by the thermocline temperature-induced pressure gradients. Obviously, the change in the temperature field is the main reason for these major differences, otherwise only the velocities along the boundaries would significantly differ.
With the exception of the tangential velocities along the lateral walls, the dynamical balance in the interior of the basin is geostrophy and the effect of momentum dissipation remains fairly small. Hence, most of the differences we can observe in the interior dynamics must be due to changes in the temperature field. In order to properly compare the velocity fields associated with the various closure schemes, we computed their diagnostic PG dynamics for the same temperature field, constructed as an average of the steady-state temperature fields of the models summarized in Table 2. The results confirm the high reliability of geostrophic balance in the interior of the basin whatever the model. However, dramatic differences occur along the boundaries: the strong horizontal tangential velocities in no-normal-flow boundary condition models have no counterpart in models with no-slip boundary conditions, while the vertical velocities are totally the opposite from one group to the other (with almost an order of magnitude difference). These astonishing differences are much stronger than in each model's steady state and proves that the boundary conditions are responsible for producing the large differences in the final steady-state temperatures.

b. Convection and vertical velocity fields
By examining the vertical velocities as a consequence of the horizontal circulation, we can infer the main differences appearing in the north-east corner: large downwelling in PGL sums up to a few Sverdrups at the base of the thermocline and contributes significantly to the MOSF intensity, but also to the heat budget of the deep water in spite of the weak vertical gradient of temperature in this region. The vertical velocity field in PGS/PGRW is amazingly regular and smooth, in comparison to the high variability in PGL, especially along the boundaries. Another major discrepancy between PGRW and PGL, which is not noticeable in the horizontal fields (although a similar consequence of the tangential velocities along the boundaries) is the upwelling in the WBC, precisely in the first grid-box along the western coast. This unfortunate feature is the well-known "Veronis effect" (Veronis 1975): strong horizontal diffusion of heat across the steeply sloping isopycnals on the west side of the Gulf Stream requires large vertical velocities (transporting heat downward) to balance it. The net upwelling along the western boundary reaches 6 Sv at 1000 m deep in PGL, in comparison to less than 2 Sv (maximum at 250 m deep) in PGRW/PGS. In PGL, this upwelling is necessary to balance the downwellings along the northern boundary and the north-east corner, in conjunction with the interior upwelling. This upwelling is also likely to contribute to the strong differences observed in the zonal overturning. The efficient reduction of these vertical transports along the boundaries is the best achievement of the no-normal-flow boundary condition, although Huang and Yang (1996) claim it is a necessary feature of frictional western boundary layers. Analyzing the results from the Community Modeling Effort, Bšning et al. (1995) consider the Veronis effect as primarily responsible for the insufficient meridional overturning rate and northward heat transport in their modeled Atlantic. Oceanic observations and evaluations of such a crucial climatic parameter are unfortunately lacking. Stommel-Arons model (1960), explaining the deep western boundary current before it was observed, is based on the planetary vorticity balance and the net upward flow at the base of the thermocline (balancing the deep-water formation). Veronis' analysis of Holland's (1971) numerical experiments points out large upwelling in the western boundary current (due to the large horizontal diffusivity) as the cause for the general downwelling at the base of the thermocline. Holland resolved this contradiction with Stommel-Arons theory in a later experiment by using a reduced horizontal diffusivity. Colin de Verdière (1989) also showed a relatively good agreement for the Sverdrup balance and Stommel-Arons's type of deep circulation in the heated region of his basin. Even with the low value of horizontal diffusivity we use here, the western upwelling remains large and strongly influences the overturning: we try to clarify the sensitivity of this numerical artifact in Huck et al. (1997).
From the high-resolution non-hydrostatic simulations of Send and Marshall (1995), a good agreement between convection regions and large-scale downwelling velocities is not necessary since convective plumes transport properties by turbulent mixing without significant net mass transports. Heat transfers by convection do not directly appear in the MOSF, although they modify the properties of the deep waters and aid in setting up east-west pressure gradients which ultimately drive the overturning. For comparison purposes, Fig. 7 presents the vertical transports at 850 m (the base of the thermocline) and the convection region across this level for PGRW and PGL. Obviously, PGRW/PGS achieve a nice agreement between downwelling and convection in the northern regions, as opposed to PGL where the noisy vertical velocity field is rarely associated with convection. We suggest that these strong vertical velocities, decoupled from static instability areas, explain a large part of the discrepancies in the MOSF and ZOSF, as well as in the bottom temperatures. Indeed the convection regions are the ones with the coldest temperatures at the surface. If the downward velocities occur in the same area, the deep waters are fed by the coldest surface waters; if strong downwellings occur elsewhere, deep water results from a mixing of water cooled by convection and water advected vertically from the surface, and then spread horizontally: finally, their temperature is higher than what can be expected from the first process only.
Comparing the deep convection regions (Fig. 4 and 7) shows an area of the northern boundary 3 times larger in PGL than PGRW, where convection is more concentrated in the western part as a consequence of the westward velocities along the northern boundaries. The north-east corner shows more convection in PGL as the large eastward current drives all the surface water there. The dynamical boundary condition along the lateral boundaries is mainly responsible for such differences in the vertical velocity fields. The simple scalings presented earlier justify an order of magnitude difference in the vertical velocities computed from the continuity equation with the no-slip boundary condition and the vorticity (or non-hydrostatic) closure. The free-slip boundary condition should not change the continuity balance much, as it affects only the normal derivative of the tangential velocity to the boundary. Although the mean vertical velocities in the interior are very similar in all the models, they only differ by a factor 3 to 5 at boundaries from no-normal-flow to no-slip boundary condition (mean vertical velocities reach more than 10-5 m s-1 along the boundaries in PGL and PG0, i.e. more than 50 times the interior average velocities). The vertical transports occurring at the boundaries are accordingly large (around 30 Sv for PGL and 40 Sv for PG0 at the depth of 2000 m, compared to 12 Sv peaking at 700 m for PGRW). Column 13 of Table 2 shows the maximum (over depth) net upwelling integrated from south to north in the one-grid box along the western boundary, hence it is a simple measure of the Veronis effect. These results suggest that the type of boundary condition makes the larger discrepancies, except in the purely geostrophic cases. We showed earlier the high correlation between this upwelling and the deep water temperature (Fig. 5d). The dramatic variability of the vertical velocity field along the boundaries generates high vertical transports of mass (and heat) that profoundly modify the zonal and meridional overturning. Since the interior dissipation is generally negligible compared to the geostrophic balance, it does not modify the vertical velocities in the interior of the mid-latitude basin. Consistent with Colin de Verdire (1988), our results confirm that the Sverdrup balance is approximately satisfied in the tropical region. In addition, the boundary condition controls the rate of upwelling occurring along the western boundary (measuring the 'Veronis effect') and these are the most intense vertical transports occurring in the basin, varying from 2 to 12 Sv. These spurious vertical velocities surely decrease the efficiency of the thermohaline circulation in two ways: by pumping cold water out of the deep WBC into the WBC ("short cut" for the deep water path), and by producing deep waters by downwelling in the north-east corner that warm the waters formed by convection in the cooler mid-northern regions. Globally, the higher mass transports, caused by this WBC-upwelling and north-east downwelling, are associated with the warmest bottom waters: this is usually correlated with lower poleward heat transport in this set of experiments. In terms of net vertical transport, the "Veronis effect" is the most important.

6. Energetics

Following Colin de Verdière (1988), we compare the energetic balance for each models' steady state. The objective is to evaluate the dissipation rate associated with the choice of momentum dissipation and boundary conditions, and the consequences on the global balance. a. The potential energy balance The potential energy equation, resulting from the basin integration of (gz) times the density equation, leads to the balance between a source term (related to the vertical mixing), two sinks (from the convection and the momentum dissipation) and a residual term (from the surface fluxes which tends to zero at steady state). These terms have the following expressions (in W m-2), with A being the total surface area of the ocean basin: where Conv(x,y,z) is the interior redistribution of density necessary to remove static instabilities (r0 a times the Convection term in the temperature equation in section 2a). Note that this value does not depend on a reference level since: The dissipation term is computed from the cross-correlation of vertical velocities and density, in order to use the same formula whatever momentum dissipation scheme is used: The residual term from the surface heat flux (depending on the reference level) is always negligible after 3000 years. We prefer to measure simply the residual surface heat flux: The kinetic energy (KE) is computed in different ways according to the numerical grid: on an A-grid, the velocities are weighted by the volume of the surrounding grid-box, but on the B-grid, either an average over the 4 corners is done to mimic the A-grid calculation, or the boundary values are weighted by one half of the volume of the grid-box. These two methods lead to significant differences which enlighten the importance of the velocities along the boundaries in the total KE. The available potential energy (APE further) is computed in two different ways: the first formula used is the one proposed by Oort et al. (1989), measuring the departure of the isopycnals from a reference state defined as the horizontal mean of the densities. A second method is also tested as it is easily implemented due to the linear equation of state and the box geometry: the reference stratification is determined by successively tracking the coldest grid-boxes in the basin and summing up their volumes to fill the basin layers, starting from the bottom. This gives a mean temperature for all the levels by an adiabatic redistribution of the ocean elements conserving their volume and temperature. Then, the APE is computed as the difference between the initial potential energy and the potential energy of the horizontally uniform reference state, using the same reference height (this arbitrary value does not influence the result since the thermal content of the two states are equal). Comparison of the two methods shows the uncertainty in the determination of the APE. Although the order of magnitude is fairly similar, the ordering of the models is not always the same, and the range of variation of APE using the second method is much smaller than with the first one (first order linearization).

b. Comparison of the energy balance for the different models
As expected from the variations in bottom water temperature, the energy diagnostics in Table 3 show large variability in the APE according to the dynamical closure. PGRW/PGS achieve the lowest values, increased by a factor of 2 to 4 (depending on the formula used) for Laplacian friction or no-slip boundary conditions, the purely geostrophic model reaches the highest APE. Understanding the ordering of APE for the different models is not trivial, as it is the opposite of the surface to bottom density contrasts which controls the reference stratification potential energy. The ordering seems correlated with the variability in the temperature field, usually dependent on the vertical velocity activity along the boundaries.
The kinetic energy (KE) interpretation is complicated by the two formulations: the B-grid value, when available, is always bigger than the A-grid one, as the averaging smoothes the largest velocities. The discrepancy between the two values is rather large when no-slip boundary conditions are used, as the vorticity closure (or non-hydrostatic with small vertical friction) generates high tangential velocities at the boundaries. Thus, a major part (60% at least) of the KE in the PGRW/PGS models comes from the boundary. When no interior dissipation is used (PG0), the highest KE is achieved. In this case only, the use of frictional vorticity closure (PG0W) to solve for lateral velocities with the no-normal-flow boundary condition introduces some dissipation and reduces the energy of the steady state. Within the same boundary condition, linear friction always induces a significantly lower KE than Laplacian viscosity. The term arising from the surface fluxes in the potential energy equation is usually less than 10-4 W m-2 after 3000 years of integration, inducing an accurate balance between the source and sink terms. The only source term for potential energy is that through vertical mixing, which we find comparable amongst each closure scheme. This is related to the similarity of surface to deep water characteristics in each model. Since the mean surface temperature needs to match the mean restoring temperature (13.5 °C) at steady-state, the relative variation of this source term is easily predicted from the variability of the mean bottom temperature (3.5±0.2 °C) over the surface to bottom difference (~10 °C), i.e. a few percent change.
The comparison of the sink terms is more interesting since various kinds of balance between convection and dissipation occur. First, PG0 shows an extreme case where the vertical mixing source is entirely dissipated by convection: purely geostrophic velocities do not produce any work and the imposed "no-slip" boundary condition satisfies a zero net dissipation (Bryan 1969; Colin de Verdire 1988). In the Laplacian cases, convection is a sink which is twice as important as momentum dissipation, which even reduces with lower viscosity. This is not changed by free-slip boundary conditions, showing that most of this dissipation happens in the interior. The linear friction model with no-slip boundary condition is very similar, with 15% more dissipation and consequently lower convection. On the opposite, PGRW and PGS loose potential energy mainly by KE dissipation (twice as important as convection). A tentative interpretation is that the vertical velocity field in these models efficiently moves the coldest water downward, having the same effect as convection (the more important rate of dissipation of energy is consistent with the lower KE with linear friction). A crude evaluation of this dissipation is possible with linear friction (as 2 eH KE) for PGR0 or PGRW (and PGS with low vertical friction), and this compares to the dissipation depending on which scheme is used to compute KE. Increasing the vertical friction coefficient in PGS makes the part of the dissipation of vertical kinetic energy more significant, reaching 10% of the horizontal one for eV=0.3. The type of boundary condition proves itself important again by its influence on the vertical velocities along the boundary. We already measured the influence of the shift from a no-slip to a no-normal-flow boundary condition with Rayleigh friction, inverting the ratio of convection and dissipation.
Comparing PG0 and PG0W is more striking, as the frictional vorticity closure for the lateral boundaries is enough to get a convection-dissipation ratio identical to PGRW. Even with no interior dissipation, the tangential velocities and associated vertical velocities along the boundaries reduce the role of convection by 60%. If our models were resolving the Rossby radius of deformation, the reduction of the convective adjustment part would acknowledge the good performance of the model in producing the right vertical velocities to resolve convection instead of parameterizing it. At coarse resolution, this argument do not allow to evaluate the models. 7. Conclusion We have compared different parameterizations of momentum dissipation and dynamical boundary conditions in the same planetary geostrophic mainframe model, with identical tracer advection, diffusion, convection and forcing. The Reynolds stress parameterizations we implemented are the traditional Laplacian friction, the higher-order biharmonic formulation, the first-order linear Rayleigh friction (in two and three dimensions) and also no momentum dissipation at all. The associated boundary conditions are either no-slip, free-slip or no-normal-flow (with a vorticity closure for the velocities tangential to the boundaries). The steady-state temperature fields and dynamics after 3000 years of integration allowed a straightforward classification of the models. The PGL group (Laplacian friction), typical of the classical Modular Ocean Model (GFDL) at coarse resolution, encompasses no-slip and free-slip Laplacian or Rayleigh friction, and free-slip biharmonic friction. This group is characterized by large vertical velocities along the boundaries, strong western boundary current (WBC) and overturning transports, a large downwelling of surface waters in the north-east corner of the basin and a large upwelling along the western boundary. Its efficiency in terms of poleward heat transport is average. The PGRW group includes the two models with no-normal-flow boundary conditions: the horizontal Rayleigh friction with vorticity closure for boundary velocities and the non-hydrostatic model with linear friction in the three directions of movement proposed by Salmon (1986, 1990). Its main features are the highest poleward heat transport, the lowest deep water temperature and a weak meridional overturning and WBC transport. The last model on its own is PG0, the purely geostrophic model with no-slip boundary conditions: it illustrates an extreme potential energy balance between vertical mixing source and convection sink, with no kinetic energy dissipation at all. Its circulation is very irregular, with a very concentrated one grid-box wide WBC and alternating northward-southward flowing bands on the east of this thermally-driven Gulf Stream. Its major deficiencies are the highest and most variable vertical velocities, the highest zonal and meridional overturnings, the warmest deep water and the lowest poleward heat transport.
These results suggest that the boundary condition has a major influence on the performance of the models, as measured by these large-scale diagnostics. Meridional and zonal overturnings are highly sensitive, almost doubling from PGRW to PG0. We found the origin of such a variability to be the vertical velocities along the boundaries: mainly, the upwelling in the northern part of the WBC (known as the 'Veronis' effect). As expected from the scaling of the continuity equation, vertical velocities along the boundaries reach 10-5 m s-1 on average with no-slip or free-slip boundary condition, producing vertical transports of up to 1 Sv per grid-box. Hence, the main effect of the no-normal-flow boundary condition is to reduce (by a factor of 3 in the mean) these huge vertical transports by allowing a horizontal recirculation of flows impinging into the coast. Otherwise, the large downwellings (in the north-east part of the basin) modify the properties of deep water formed by convection while the large upwelling in the WBC provides a short-cut to the thermohaline loop, strongly reducing the intensity of the deep WBC as it travels southward and hence, the total poleward heat transport. Since the interior of the basin satisfies a geostrophic balance, weak vertical transports occur there and the role of the boundaries is fundamental in the global exchange of mass. The sensitivity to the interior momentum dissipation (mixing by Laplacian or biharmonic, sink of energy by linear friction) is of second-order influence compared to no-dissipation influence. The temperature fields vary according to these strong dynamical discrepancies along the boundaries. Since the mean surface temperature is fixed by the restoring boundary condition, we concentrated on the deep water characteristics and observed differences of up to 0.5 °C, depending on the closure scheme: this is easily comparable to the influence of tracer mixing coefficients. These bottom temperatures have a strong influence on the ability of the model to transport heat poleward, in negative correlation with the meridional overturning and the deep water temperature. Since the increase in meridional overturning and the warming of deep water is highly correlated to the Veronis effect, we confirm recent studies enlightening the short-cut of the thermohaline loop as primarily responsible for the weak North Atlantic poleward heat transport in global models (Bšning et al. 1995) although the real ocean upwelling in the Gulf Stream is not known. The role of dynamical boundary conditions in coarse-resolution thermohaline circulation models is more important than the interior momentum dissipation. Laplacian friction produces a more satisfying interior circulation, in better agreement with geostrophy, but generates huge vertical transports along the lateral boundaries (especially upwelling in the western boundary current and downwelling in the north-east corner). Rayleigh friction, with a vorticity closure whose primary effect is to reduce vertical velocities along the boundaries by allowing horizontal recirculation, induces a more efficient thermohaline circulation (colder deep water, weaker meridional overturning but higher poleward heat transport). This parameterization lacks physical justifications and is not as satisfying as the Laplacian in terms of the interior geostrophic and Sverdrup balance. It is however an interesting alternative to implement no-normal-flow boundary conditions, since free-slip and no-slip boundary conditions are shown here to lead to very similar circulations.
Hence, the closure scheme is an important and overlooked aspect of climate modeling studies, comparable to tracer mixing parameterization in its global effects. A proper closure for the dynamical boundary conditions based on the physics of the real ocean near lateral boundaries may dramatically improve the results of our OGCMs. However, wind-forcing and a realistic topography may reduce the influence of boundary conditions by constraining more strongly the deep-water formation regions.

Acknowledgments. This research was funded through a scholarship awarded to TH from the French Ministry of Research and Space and operating grants awarded to AJW from NSERC, AES, CICS and the NOAA Scripps-Lamont Consortium on the Ocean's Role in Climate. We would like to thank A. Fanning for many comments and fruitful discussions about this work, as well as an anonymous reviewer for numerous thoughtful suggestions.

REFERENCES

Boening, C. W., W. R. Holland, F. O. Bryan, G. Danabasoglu, and J. C. McWilliams, 1995: An overlooked problem in model simulations of the thermohaline circulation and heat transport in the Atlantic ocean. J. Climate, 8, 515-523. Bryan, F., 1987: Parameter sensitivity of primitive equation ocean general circulation models. J. Phys. Oceanogr., 17, 970-985. Bryan, K., 1962: Measurements of meridional heat transport by ocean currents. J. Geophys. Res., 67, 3403-3414. Bryan, K., 1969: A numerical method for the study of the circulation of the world ocean. Journal of computational physics, 4, 347-376. Bryan, K., 1991: Poleward heat transport in the ocean - A review of a hierarchy of models of increasing resolution. Tellus, 43 AB, 104-115. Colin de Verdière, A., 1986: On mean flow instabilities within planetary geostrophic equations. J. Phys. Oceanogr., 16, 1981-1984. Colin de Verdière, A., 1988: Buoyancy driven planetary flows. J. Mar. Res. , 46, 215-265. Colin de Verdière, A., 1989: On the interaction of wind and buoyancy driven gyres. J. Mar. Res., 47, 595-633. Danabasoglu, G., and J. C. McWilliams, 1995: Sensitivity of the global ocean circulation to parameterizations of mesoscale tracer transports. J. Climate, 8, 12, 1967-2987. Danabasoglu, G., J. C. McWillliams, and P. R. Gent, 1994: The role of mesoscale tracer transports in the global ocean circulation. Science, 264, 1123-1126. Gent, P. R., and J. C. McWilliams, 1990: Isopycnal mixing in ocean circulation models. J. Phys. Oceanogr., 20, 150-155. Gent, P. R., and J. C. McWilliams, 1996: Eliassen-Palm fluxes and the momentum equation in non-eddy-resolving ocean circulation models. J. Phys. Oceanogr., 26, 2539-2546. Greatbatch, R. J., and K. G. Lamb, 1990: On parameterizing vertical mixing of momentum in non-eddy resolving ocean models. J. Phys. Oceanogr., 20, 1634-1637. Haney, R., 1971: Surface thermal boundary condition for ocean circulation models. J. Phys. Oceanogr., 1, 241-248. Hasselmann, K., 1982: An ocean model for climate variability studies. Prog. Oceanog., 11, 69-92. Hirst, A. C., and T. J. McDougall, 1996: Deep-water properties and surface buoyancy flux as simulated by a Z-coordinate model including eddy-induced advection. J. Phys. Oceanogr., 26, 1320-1343. Holland, W., 1971: Ocean tracer distributions - a preliminary numerical experiment. Tellus, 23, 371-392. Huang, R. X., and J. Yang, 1996: Deep-wtaer upwelling in the frictional western boundary layer. J. Phys. Oceanogr., 26, 2243-2250. Huck, T., Colin de Verdire, A., and A. J. Weaver, 1997: The effect of momentum dissipation parameterizations in OGCM: Geostrophy, Sverdrup balance and Veronis effect. in preparation. Huck, T., 1997: Modeling the large-scale ocean thermohaline circulation: Analysis of its decadal variability. Ph. D. thesis manuscript, University of Brest, France. Killworth, P., 1985: A two-level wind and buoyancy driven thermocline model. J. Phys. Oceanogr., 15, 1414-1432. McCartney, M. S., and L. D. Talley, 1982: The subpolar mode water of the North Atlantic ocean. J. Phys. Oceanogr., 12, 1169-1188. Marotzke, J., 1991: Influence of convective adjustment on the stability of the thermohaline circulation. J. Phys. Oceanogr., 21, 903-907. Munk, W. H., 1950: On the wind-driven ocean circulation. J. Meteorology, 7, 2, 79-93. Oort, A. H., S. C. Ashcher, S. Levitus, and J. P. Peixoto, 1989: New estimates of the available potential energy in the world ocean. J. Geophys. Res., 94, 3187-3200. Pacanowski, R. C., 1995: MOM 2 Documentation, User's Guide and Reference Manual. GFDL Ocean Technical Report #3, 232 pp. Pacanowski, R., Dixon, K., and A. Rosati, 1991: The GFDL Modular Ocean Model, Users Guide Version 1.0. GFDL Ocean Group Technical Report #2, 46 pp. Phillips, N. A., 1963: Geostrophic motion. Rev. Geophysics Space Physics, 1, 123-176. Rahmstorf, S., 1994: Multiple equilibria of an OGCM caused by different convection patterns. Ocean Modelling, 103 (unpublished manuscript). Rahmstorf, S., 1993: A fast and complete convection scheme for ocean models. Ocean Modelling, 101 (unpublished manuscript). Robitaille, D. Y., and A. J. Weaver, 1995: Validation of sub-grid scale mixing schemes using CFCs in a global ocean model. Geophys. Res. Letters, 22, 2917-2920. Salmon, R., 1986: A simplified linear ocean circulation theory. J. Mar. Res., 44, 695-711. Salmon, R., 1990: The thermocline as an "internal boundary layer". J. Mar. Res., 48, 437-469. Samelson, R. M., and G. K. Vallis, 1997a: A simple friction and diffusion scheme for planetary geostrophic basin models. J. Phys. Oceanogr., 27, 186-194. Samelson, R. M., and G. K. Vallis, 1997b: Large-scale circulation with small diapycnal diffusion: The two-thermocline limit. J. Mar. Res., 55, 223-275. Semtner, A. J., 1986: Finite-difference formulation of a world ocean model. Proceedings of the NATO Advanced Study Institute on Advanced Physical Oceanographic Numerical Modelling, D. Reidel Publishing Co, Dordrecht Semtner, A. J., and R. M. Chervin, 1988: A simulation of the global ocean circulation with resolved eddies. J. Geophys. Res., 93, 15,502-522. Send, U., and J. Marshall, 1995: Integral effects of deep convection. J. Phys. Oceanogr., 25, 855-872. Stewart, R. W., 1989: The no-slip constraint and ocean models. Atmos.-Ocean, 27, 542-552. Stommel, H., and A. B. Arons, 1960: On the abyssal circulation of the world ocean. II: An idealized model of the circulation pattern and amplitude in oceanic basins. Deep-Sea Res., 6, 217-233. Stommel, H., 1948: The westward intensification of wind-driven ocean currents. Transactions, American Geophysical Union, 29, 2, 202-206. Straub, D. N., 1996: An inconsistency between two classical models of the ocean's buoyancy driven circulation. Tellus, 48A, 477-481. Veronis, G., 1973: Large scale ocean circulation. Adv. Applied Mechanics, 13, 1-92. Veronis, G., 1973: Models of world ocean circulation: 1. Wind-driven, two layer. J. Mar. Res., 31, 228-289. Veronis, G., 1975: The role of models in tracer studies. Numerical models of the ocean circulation, 133-146. Nat. Acad. of Sciences, Washington DC 14pp. Weaver, A. J., and M. Eby, 1997: On the numerical implementation of advection schemes for use in conjunction with various mixing parameterizations in the GFDL Ocean Model. J. Phys. Oceanogr., 27, 369-377. Weaver, A. J., and E. S. Sarachik, 1991: Evidence for decadal variability in an ocean general circulation model: an advective mechanism. Atmos.-Ocean, 29, 197-231. Winton, M., and E. S. Sarachik, 1993: Thermohaline oscillations induced by strong steady salinity forcing of ocean general circulation models. J. Phys. Oceanogr., 23, 1389-1410. Winton, M., 1993: Numerical investigations of steady and oscillating thermohaline circulations. Ph. D. thesis, University of Washington, 155 pp. Zhang, S., C. A. Lin, and R. J. Greatbatch, 1992: A thermocline model for ocean-climate studies. J. Mar. Res., 50, 99-124.

FIGURE CAPTIONS

FIGURE 1. Minimum, mean bottom and mean basin temperature (°C) for PGRW (solid lines) and PGL (dashed lines) as a function of horizontal tracer diffusivity (in a log-log plot). For each model, the upper curve is the mean basin temperature, the middle one is the mean bottom temperature (at 4225 m deep) and the lower curve is the minimum temperature, usually achieved from the surface to the bottom by convection. Note the shift between the 2 lowest curves for PGL compared to PGRW (related to the variability in the deep water temperatures) and the inversion at low diffusivity due to numerical problems (high Peclet number).
FIGURE 2. Zonally-averaged temperature in the thermocline (°C) and maximum depth of convection (dashed line) for (a) PGRW, (b) PGL and (c) PG0. The thermocline becomes more diffuse, mainly because of the warmer deep water temperatures.
FIGURE 3. Meridional overturning streamfunction (Sv) and maximum depth of convection (dashed line) for (a) PGRW, (b) PGL and (c) PG0. Note the large variation in the maximum intensity and its position, as well as the more localized convection regions in PGRW.
FIGURE 4. Zonal overturning streamfunction (Sv), meridionally-averaged isotherms (dashed lines) and maximum convection depth over latitude (shaded area) for (a) PGRW, (b) PGL and (c) PG0. Note the slope of the isotherms in the western boundary current for PGL and the narrow region of deep convection in PGRW.
FIGURE 5. Cross correlations of diagnostics for models 3 to 16 in Table 2: maximum advective poleward heat transport function of meridional overturning (a) and mean bottom temperature (b); meridional overturning (c) and western upwelling (d) function of mean bottom temperature. r is the correlation coefficient (the 95% confidence level for the 14 models used is 0.52) and the solid line is the linear regression solution.
FIGURE 6. Mean temperature (°C) and horizontal velocities in the upper 850 m (thermocline) for (a) PGRW, (b) PGL and (c) PG0. The scale for velocities is 2 cm s-1 per grid spacing. Note the western boundary current overshoot, characteristic of the Laplacian viscosity.
FIGURE 7. Vertical transports (Sv) at the base of the thermocline (850 m) and regions of convection below 250 and 850 m (dashed lines) for (a) PGRW, (b) PGL and (c) PG0. The downwelling regions are shaded and the appropriate vertical transports are then downward. Note the large and irregular transports along the boundaries in PGL and PG0 compared to PGRW.

TABLE CAPTIONS

TABLE 1. Models summary: momentum dissipation and boundary conditions. All dynamical boundary conditions imply no mass transport across the boundaries: in addition, no-slip adds the constraint of zero tangential velocities along the boundaries, while free-slip requires the normal derivative of tangential velocities to vanish along the boundary (no stress); no-normal-flow allows tangential velocities along the boundaries, that need to be determined through an additional equation.
TABLE 2. Steady-state diagnostics for the different numerical experiments. Rows 1 and 2 measure the influence of the vertical viscosity with the implicit vertical mixing option for convection in the MOM. Rows 3 to 16 implement Rahmstorf's full convection scheme and no vertical viscosity. Column 1 is the experiment number (used in Fig. 5). Column 2 is the model name. Column 3 is the mean basin temperature, column 4 the averaged bottom temperature from 3950 to 4500 m, column 5 the minimum temperature achieved at the surface and extending to the lower levels by deep convection. Column 6 is the maximum meridional overturning (as defined in section 4b). Column 8 (7) is the maximum (advective) poleward heat transport (section 4c). Column 9 is the maximum (thermally driven) Gulf Stream transport, and the latitude of the maximum is given in column 10. Columns 11 and 12 are the minimum and maximum of the zonal overturning streamfunction (defined in section 4b). Column 13 is the maximum (over depth) of the vertical transport in the one-grid-box wide region extending along the whole western boundary. Finally, column 14 is the non-dimensional mean kinetic energy in steady-state.
TABLE 3. Energy balance for different steady-states. Column 1 is the model name. Columns 2 and 3 give the kinetic energy computed on an A and B grid Columns 4 and 5 give estimations of the available potential energy with the Oort formula or with the method given in section 6a. Columns 6, 7 and 8 are the potential energy balance with the vertical mixing source term and the convective and dissipative sinks. The numbers in parenthesis refer to the estimation of vertical or horizontal linear dissipation of kinetic energy computed by 2 eV KEVertical and 2 eH KE. Finally, the residual mean surface heat flux is given in column 9. Details of the computations are given in section 6a.

Table 1. Models summary: momentum dissipation and boundary conditions. Momentum Lateral Horizontal Number Dissipation Boundary Grid of Horizontal Vertical Conditions Type levels MOM 0 no slip B 15 + non-linear terms PGL 0 no slip B 15 PGLA 0 no slip A 15 PGLslip 0 free slip B 15 PGBslip 0 free slip B 15 PG0 0 0 no slip B 15 PG0slip 0 0 free slip B 15 PG0W 0 0 no normal flow B 15 tangential velocities via frictional vorticity closure PGR0 0 no slip B 15 PGRslip 0 free slip B 15 PGRW 0 no normal flow B 15 tangential velocities via frictional vorticity closure PGRW90 0 no normal flow B 90 tangential velocities via frictional vorticity closure PGR4 0 no normal flow B 15 biharmonic dissipation is added to the tracer equation PGS no normal flow B 90

Table 2. Steady-state diagnostics for the different numerical experiments. exp ## Model Viscosity Mean Temp. (¡C) Bottom Temp. (¡C) Min. Temp. (¡C) max. MOSF (Sv) max. a PHT (PW) max. PHT (PW) max. G S (Sv) pos. max. (¡N) min. ZOSF (Sv) max. ZOSF (Sv) West. Upw. (sv) mean KE (.1 J m-3) 1 MOM Av=10-3 4.228 3.531 3.461 15.42 .2204 .2308 15.0 42 -4.08 7.62 5.06 .2562 2 MOM Av=0 4.228 3.531 3.462 15.41 .2204 .2308 15.0 42 -4.07 7.63 5.05 .2562 3 MOM fullconv 4.230 3.533 3.464 15.39 .2204 .2308 15.0 43 -4.04 7.87 5.05 .2600 4 PGL 4.237 3.539 3.470 15.45 .2219 .2321 15.3 42 -3.84 7.99 5.14 .2752 5 PGLA 4.251 3.547 3.478 15.47 .2203 .2304 14.9 43 -3.68 8.12 5.84 .2490 6 PGLslip 4.189 3.477 3.405 11.14 .2258 .2360 16.0 38 -3.74 8.64 4.90 .3350 7 PGBslip 4.238 3.509 3.434 10.70 .2235 .2341 15.0 40 -3.54 9.38 5.87 .4716 8 PG0 4.530 3.748 3.653 15.57 .2133 .2248 19.2 33 -2.60 12.10 12.30 .9690 9 PG0slip 4.262 3.529 3.449 10.91 .2238 .2346 16.6 34 -3.56 9.67 6.30 .5733 10 PG0W 4.127 3.413 3.388 7.57 .2477 .2602 14.3 30 -2.03 3.70 2.95 .6940 11 PGR0 4.171 3.513 3.441 14.25 .2142 .2236 11.8 42 -3.46 6.89 4.75 .1930 12 PGRslip 4.112 3.458 3.402 12.49 .2175 .2268 12.0 42 -3.07 7.12 3.78 .1880 13 PGRW 4.056 3.404 3.389 9.00 .2374 .2471 10.5 40 -2.30 2.79 1.94 .3590 14 PGR4 4.200 3.602 3.586 7.98 .2287 .2389 9.71 37 -0.44 2.14 0.63 .2968 15 PGRW90 4.083 3.416 3.401 9.09 .2378 .2489 10.4 40 -2.21 2.85 1.93 .3546 16 PGS ev=0.03 4.080 3.406 3.392 8.71 .2375 .2486 10.3 39 -1.99 2.75 1.77 .3360 17 PGS ev=0.3 4.061 3.346 3.329 7.63 .2330 .2442 9.5 38 -1.36 2.45 1.46 .2909

Table 3. Energy balance for different steady-states. Model KE A-grid (J/m2) KE B-grid (J/m2) APE Oort (105 J/m2) APE Ref. (105 J/m2) DPE Vert. Mixing (10-3 W/m2) -DPE Convection (10-3 W/m2) -DPE Dissipation (10-3 W/m2) Residual SHF (W/m2) MOM fullconv 104 117 8.75 3.51 1.99 1.33 0.66 4.84´10-5 PGL 110 124 7.96 3.47 1.99 1.34 0.65 1.70´10-5 PGLA 112 7.94 3.51 2.00 1.37 0.63 1.20´10-5 PGLslip 144 150 5.59 3.13 2.00 1.35 0.65 3.69´10-5 PGBslip 187 212 8.15 3.21 1.96 1.39 0.56 5.70´10-3 PG0 212 436 13.00 4.23 1.91 1.91 0.00 1.98´10-6 PG0slip 203 258 9.35 3.28 1.95 1.43 0.53 1.85´10-5 PG0W 217 312 2.26 1.96 2.02 0.72 1.30 4.55´10-4 PGR0 73 87 7.82 3.28 2.00 1.23 0.77 (0.76) 2.10´10-5 PGRslip 82 85 6.11 2.94 1.97 1.14 0.83 (0.75) 4.47´10-5 PGRW 108 162 2.02 1.72 2.02 0.64 1.38 (1.42) 3.81´10-4 PGR4 116 134 2.076 1.80 1.94 0.82 1.12 (1.17) 2.00´10-3 PGS ev=0.03 104 151 2.008 1.91 1.98 0.65 (0.07) 1.33 (1.32) 3.13´10-4 PGS ev=0.3 96 131 2.014 1.99 1.99 0.68 (0.23) 1.32 (1.15) 5.15´10-4