Decadal variability of the thermohaline circulation in ocean models.

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

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


submitted to Journal of Physical Oceanography: June 25, 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. THE MODELS

3. PARAMETER SENSITIVITY ANALYSIS
a. Methodology
b. Spherical vs. Cartesian coordinates
c. The weak influence of beta
d. Influence of the Coriolis parameter on Cartesian f-planes
e. The damping effect of the horizontal diffusion
f. The driving effect of the vertical diffusion
g. The weak influence of the convective adjustment
h. Influence of the momentum dissipation parameterization and associated boundary conditions
i. The damping influence of friction
j. Influence of the horizontal resolution
k. The sensitivity to the vertical discretization
l. Influence of the forcing amplitude
m. Influence of the basin width

4. THE ROLE OF THE BOUNDARIES
a. Motivation
b. The eastern boundary
c. The polar boundary
d. The western boundary

5. HINTS FOR THE MECHANISM IN THE TRANSITION FROM STEADY TO OSCILLATORY STATES
a. Shift in the mean state when restoring boundary conditions are changed to fixed fluxes
b. Change in the variability patterns as a function of the forcing field
c. The regions with maximum anomaly growth rate
d. A geostrophic amplification of the western anomalies
e. The meridional overturning adjustment to the density anomalies
f. A simple box-model analogy for the oscillation

6. ROBUSTNESS OF THE VARIABILITY TO VARIOUS PARAMETERIZATIONS AND FORCING

7. CONCLUSION

APPENDIX: GEOSTROPHIC SCALING OF OVERTURNING STREAMFUNCTION UNDER CONSTANT FLUX BOUNDARY CONDITION

REFERENCES

TABLE CAPTIONS

FIGURE CAPTIONS



ABSTRACT

---

Intrinsic modes of decadal variability driven by constant surface buoyancy fluxes are analyzed using a box-geometry ocean model. A complete parameter sensitivity analysis of the oscillatory behavior is carried out with respect to the spherical/Cartesian geometry, the beta-effect, the Coriolis parameter, the parameterization of momentum dissipation and associated boundary conditions and viscosities, the vertical and horizontal diffusivities, the convective adjustment parameterization, the horizontal and vertical model resolution, the forcing amplitude and the basin width. The oscillation stands out as a robust geostrophic feature whose amplitude is mainly controlled by the horizontal diffusivity. Attempts to reproduce the variability in zonally-averaged 2D models being unsuccessful, the 3D adjustment processes seem necessary. Since the b-effect has no major influence, it is not the classical Rossby waves that play the leading role in the mechanism (Winton 1996). Then, various experiments with different geometry and forcing are conducted to test the importance of numerical boundary waves in sustaining the oscillation. The results suggest that only the western boundary is crucial (Greatbatch and Peterson 1996). The analysis of the variability patterns differentiates two types of oscillatory behavior: temperature anomalies propagating geostrophically westward in the eastward jet (northern part of the basin) and inducing an opposite anomaly in their wake; temperature anomalies in the north-west corner which respond to the western boundary current transport changes, but reinforce geostrophically this change and build the opposite temperature anomaly in the east, which finally reverse the meridional overturning anomaly (and thus the anomalous western boundary current transport). The analysis of the transition from steady to oscillatory states, as well as the comparison of the variability under various forcing fields, suggest that the variability is triggered in the regions of strongest surface cooling. Finally, we propose a simple box-model analogy that captures the crucial phase-shift between meridional overturning and north-south density gradient anomalies at these decadal time scales.

---

1. Introduction

Analysis of hydrographic and other data has recently revealed that the coupled ocean-atmosphere system is variable on decadal time-scales (Deser and Blackmon 1993; Kushnir 1994; Hansen and Bezdek 1996; Reverdin et al. 1997). This has important ramifications for global climate change since before one can unambiguously detect a signal of anthropogenic global warming, one must acquire an understanding of the free modes of variability of the climatic system. Furthermore, the analysis of the present state of the coupled system in terms of these modes reveals the potential for some predictability on a few years (Griffies and Bryan 1997). Coupled atmosphere-ocean models appear to be able to reproduce observed patterns of interdecadal variability involving the thermohaline circulation (Delworth et al. 1993). Latif and Barnett (1996) propose a coupled mechanism for interdecadal climate variability involving the slow adjustment of the ocean circulation to changes in the wind patterns induced by sea surface temperature anomalies, whereas Wohlleben and Weaver (1995) and Timmermann et al. (1996) outline a mechanism involving the coupling of the thermohaline circulation with atmospheric stormtracks and freshwater export out of the Arctic. Frankignoul et al. (1996) suggest that a white stochastic wind forcing can produce dominant decadal periods of ocean variability, related to the time long baroclinic Rossby waves propagate across the basin. On the other hand, Weaver and Sarachik (1991) and Yin and Sarachik (1995) propose an advective-convective mechanism for the decadal oscillations found in idealized ocean models forced with mixed boundary conditions, involving the circulation of salinity anomalies in the sub-polar gyre. Numerous ocean models forced by fixed flux boundary conditions generate purely oceanic decadal-scale variability (Huang and Chou 1994-hereafter HC94; Chen and Ghil 1995; Greatbatch and Zhang 1995; Winton 1996-hereafter W96; Greatbatch and Peterson 1996-hereafter GP96). Even a zonally-averaged ocean model coupled to a 2-level atmospheric model exhibits variability on the interdecadal time scale (Saravanan and McWilliams 1995). However, W96 reports no oscillation in a zonally-averaged ocean model forced by constant surface heat fluxes, although the same fluxes induce decadal oscillations when applied to a fully three dimensional model.
W96 reproduced the decadal oscillation found in ocean general circulation models on an f-plane, thus discarding classical Rossby wave (based on the b-effect) from the mechanism. Nevertheless, we will show that the variability patterns are quite different on f- and beta- planes. GP96, in agreement with W96, propose an explanation based on the propagation of frictional boundary waves (the viscous analogue of Kelvin waves when no time-derivative is retained in the momentum equation). They suggested that these waves are sufficiently slowed along the weakly stratified polar boundaries, where convection takes place, to give rise to decadal periods. We will analyze this mechanism and suggest that it fails in some cases and so is not a general result.
The relatively long time-scale of decadal oscillations (compared to atmospheric thermal inertial time scales) suggests a necessary oceanic contribution. Our aim is therefore to study the non-coupled oceanic variability in order to simplify the understanding of its oscillatory behavior. We concentrate here on the purely thermally-forced ocean circulation for simplicity and in order to get the essence of the mechanism for decadal variability in ocean models.
Most of the numerical experiments we perform are based on the planetary geostrophic equations, but direct comparisons with the Geophysical Fluid Dynamic Laboratory Modular Ocean Model (GFDL MOM) validate our choice of equations for the relatively coarse resolution used. For each experiment, we use the same geometry, parameters and forcing but allow for different parameterizations of momentum dissipation and associated boundary conditions in order to discard any conclusions that would be specific to one model. We show that decadal oscillations are easily found under zonally-uniform constant surface heat flux (varying with latitude) but also with more realistic westerly intensified fluxes diagnosed at the equilibrium under restoring boundary conditions. In order to isolate the mechanism driving the variability, we carry out a complete parameter sensitivity study, similar to the work of HC94 for the haline circulation. After discarding convection and the b-effect from the mechanism, we show that the oscillation appears to be primarily a geostrophic perturbation, very sensitive to the damping by horizontal diffusivity, and more likely to be generated by a strong circulation. The role of the boundaries is then analyzed: in comparison with GP96, we extend the ocean basins with buffer regions where no forcing is applied. Similar conclusions apply: only the western boundary seems necessary to sustain the oscillations. However, additional experiments with forcing symmetric in latitude (that removes the weakly stratified northern boundary) suggests that a boundary wave propagation mechanism is not always appropriate. Finally, we try to determine the origin of the variability in the analysis of the transient phase from a steady circulation to a regular oscillatory behavior, when restoring boundary conditions are changed to diagnosed fixed fluxes. The source of the variability appears to be related to the regions of strongest surface cooling. The amplification of such anomalies in the vicinity of the western boundary appears geostrophically driven, but we associate the oscillation period to the whole meridional overturning adjustment lag to density anomalies. Finally, we propose a simple box model formulation that captures the phase-shift between the strength of the meridional cell and the north-south density gradient anomalies.
The rest of the paper is organized as follows: section 2 briefly describes the models used, section 3 reports for the parameter sensitivity analysis while section 4 looks at the role of the boundaries. Section 5 provides hints to the variability mechanism through the transition from restoring to fixed flux boundary conditions. The robustness of the variability to various parameterizations and processes is discussed in section 6 and conclusions are given in section 7.

2. The models

The models we use are described in detail in Huck et al. (1996-H96 hereafter). They are based on the Cartesian coordinates planetary geostrophic equations (Phillips 1963; Hasselmann 1982; Colin de Verdi�re 1988, 1989; Zhang et al. 1992; Winton 1993) and rely on the hydrostatic and Boussinesq approximations. A version of these models was also developed in spherical coordinates to assess the validity of our Cartesian coordinates assumption. The GFDL MOM (Pacanowski 1995), based on the primitive equations, was also used to validate the planetary geostrophic assumption (see H96 for details). The geometry has not been modified from H96: a flat-bottomed b-plane centered at 40�N, extending from 20�N to 60�N, 5120 km wide and 4500 m deep. The resolution is similar in most of the cases: 160km�160km horizontally, 15 levels vertically (respectively 50�3, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550�3 m deep). Wind-stress is neglected in order to simplify the analysis of the purely thermohaline variability. Traditional sub-grid-scale Laplacian parameterizations are used for tracer diffusion (700 and 10-4 m2/s respectively for the horizontal and vertical diffusivities) along with the convective adjustment scheme of Rahmstorf (1993). For momentum dissipation, no vertical mixing is used in most of the cases, since the high horizontal friction needed to resolve the frictional boundary layers within our coarse grid has a much greater influence (see section 3i). Horizontal viscosity is parameterized either by the usual Laplacian operator (model PGL, with a 1.5�105 m2/s viscosity) or by a simpler linear drag (models PGR0 and PGRW, with a friction coefficient set to 6% of the rotation rate of the earth, i.e. 4.4e-6 s-1).
Fixed flux (Neumann) surface boundary conditions are used, with the "linear" forcing of section 3 consisting of zonally-uniform surface heat fluxes varying linearly with latitude from 45 W m-2 at 20�N to -45 W m-2 at 60�N. One might argue that using temperature as the only density variable with a flux boundary condition does not match with the rather fast relaxation of SST anomalies under real atmospheric forcing. Nevertheless, it has been shown that restoring boundary conditions produce unrealistic feedbacks (Zhang et al. 1993; Power and Kleeman 1994; Rahmstorf and Willebrand 1995). Actually, the large spatial scale of the SST anomalies we observe here will have a significant influence on the atmospheric temperature and thus, dramatically lengthen the typical radiative damping time scale (Seager et al. 1995). Furthermore, Greatbatch and Zhang (1995) note that in coupled models, changes in the poleward heat transport determine predominantly the local heat content in the convection regions. Here, temperature is utilized as a simple and meaningful way to deal with density and the heat forcing should be considered as a global buoyancy one. The mixed boundary conditions traditionally used in large-scale ocean modeling are split into the thermal relaxation that drives principally the thermohaline cell, and fixed freshwater flux, allowing more freedom for salinity anomalies to develop and play a leading role regarding the variability (or catastrophe) of the circulation. Thus, analysis of decadal oscillations within mixed boundary conditions have focused on the role of salinity anomalies (Weaver and Sarachik 1991; Lenderink and Haarsma 1996).
Most of the models we use implement a no-slip boundary condition: PG0 for Purely Geostrophic and PGR0 for Planetary Geostrophic with Rayleigh friction (Zhang et al. 1992), PGL for Planetary Geostrophic with Laplacian friction (Colin de Verdi�re 1988; Winton 1993). Two models depart fundamentally from this group: implementing a no-normal-flow boundary condition, they rely on a vorticity balance for solving the tangential velocities allowed along the boundaries. Since the frictional vorticity equation incorporates the continuity equation through the stretching term (f dw/dz), this closure efficiently reduces the vertical velocities along the lateral boundaries and enhances the horizontal recirculation of flows impinging into coasts, profoundly influencing the whole overturning and deep water properties. Curiously, the associated models PGRW-for Planetary Geostrophic with Rayleigh friction and Winton's vorticity closure (Winton 1993; Winton and Sarachik 1993)- and PGS-for Planetary Geostrophic with a 3-dimensional linear friction relaxing the hydrostatic approximation (Salmon 1986, 1990)- are the most reluctant to sustain variability, perhaps because of the weak overturning they produce. As reported by Greatbatch and Zhang (1995) and HC94, respectively, for the thermally and salinity driven circulations, zonally uniform buoyancy fluxes (varying with latitude) can induce decadal oscillations. Cai et al. (1995) added that a small zonal redistribution of surface buoyancy fluxes diagnosed at equilibrium under restoring boundary condition triggers interdecadal variability. In fact, in most of our models we observed that zonal redistribution is not necessary (GP96). For some of the experiments described in H96, the shift from restoring boundary conditions to fixed fluxes diagnosed at the equilibrium did not change the equilibrium. However, the same fixed-flux forcing applied to an ocean at rest (uniform temperature) produces perpetual oscillations in most of the models. This clearly suggests that initial conditions are important. Furthermore, a minimal temperature perturbation to the equilibrium state (0.1�C over 500 km � 500 km and 100 m deep) triggers oscillations under the same diagnosed fixed fluxes in all the models except PGRW. Obviously, there is a subtle combination of parameters and forcing that gives rise to different model results (Cai et al. 1995; GP96). Many of these previous studies experiments under fixed flux were only initialized with the exact equilibrium state. The question whether the oscillation behavior depends on the initial state (same model, parameters and forcing) is a fundamental issue: Numerous experiments we carried out suggest that if the model does not remain in equilibrium upon switching to fixed flux boundary conditions (that is rather unlikely in the real ocean, considering the atmospheric synoptic activity or the radiative seasonal cycle), then the oscillations that get fully developed after a few hundred years have the same characteristics for the same parameters and forcing, whatever the initial temperature field was (like a mode of the mean circulation). The statistical characteristics of the oscillations that we discuss in the following are therefore robust, regardless of the initial state (see for instance exp. 75 vs. 76).
The variability occurring in these models under constant surface heat fluxes (or under restoring boundary conditions with a long relaxation time-scale, see section 6) is concentrated in the thermocline (upper 400 m), where temperature anomalies reach a few degrees. The pattern of variability is rather complex and varies widely depending on the forcing (section 5b), but remains usually concentrated in the north-west corner. In the deeper layers, the variability is greatly reduced and concentrated in the convective regions along the polar boundary. The consequences of these temperature deviations on the model dynamics is always important in terms of meridional overturning (many Sverdrups) or poleward heat transport (0.05 to 0.2 petawatt). A typical evolution of the upper 250 m temperatures and velocities over a period is discussed in section 5d; the oscillation periods produced by the different models are in the interdecadal range (15-30 years).

3. Parameter sensitivity analysis

a. Methodology
We begin by comparing the variability found in models using either spherical or Cartesian coordinates. Since no major differences arise from the use of either coordinate system, we use Cartesian coordinates in all subsequent experiments in order to study the influence of b, the Coriolis parameter (on f-planes), the horizontal and vertical diffusion, the convective adjustment, the momentum dissipation, the horizontal and vertical discretization, the basin width and the forcing amplitude. The forcing for all these experiments is the "linear" surface heat flux described in the last section, where our control run was discussed. The numerous experiments are summarized in Table 1, while Table 2 gives various statistics for the mean circulation and its variability. The statistics are computed once the oscillatory state is regular, i.e. when the diagnostic values averaged over an oscillation period remain constant (usually after 2000 years). The temperature field is initialized at 4�C and the zero area-average of the surface heat fluxes implies that the mean basin temperature does not drift in time. HC94 carried out a similar parameter sensitivity study for the haline circulation under the natural boundary condition for salt proposed by Huang (1993), which drives a very intense overturning cell sinking at low latitudes. Here, we complete this work in the context of the thermally driven circulation, with traditional surface fluxes inducing a direct overturning cell of realistic amplitude. Therefore our conclusions may differ and the comparison with their results is done when the same parameter is considered.
Different schemes for momentum dissipation and associated boundary conditions are used (H96) to test the robustness of the conclusions we draw about the effect of processes and parameters. Since the oscillation seems very sensitive to the interaction of model parameterizations, resolution, parameters and forcing, great care has been taken when different models produced different behavior: only general conclusions are drawn. For instance, explicit horizontal diffusion was chosen in relation with centered difference advection scheme to control numerical diffusion (as opposed to upstream scheme where diffusion is proportional to velocity and resolution). In most cases, the oscillatory behavior is robust in a large range of parameters and has a significant amplitude in terms of SST and heat transport anomalies. To quantify our model inter-comparison, we define oscillation periods and indexes that do not depend on the various dynamics. When the oscillation is regular and quasi-periodic, measuring the period is straightforward. However, some cases are not regular: the estimated period is then computed via Fourier analysis of the temperature time series (a '~' is then printed before the period value in Table 2): In all the cases an accurate maximum in the energy-frequency spectrum determines the most energetic oscillation period. The index needs to be defined on the temperature field, the only common ground between models. Therefore, we chose the basin average of the standard deviation of temperature over one period, if the oscillation is regular, otherwise over a long time compared to the estimated period:
Since this index is difficult to interpret regarding its consequences on the circulation (depending on the various dynamics), traditional diagnostics and their standard deviations are added in Table 2.

b. Spherical vs. Cartesian coordinates
Spherical geometry is necessary for realistic modeling of the oceanic basin-scale thermohaline circulation. However, for the purpose of understanding, Cartesian coordinates are far more tractable and will be used for the further parameter sensitivity analysis. Thus, we first justify the validity of Cartesian coordinates for studying these oscillations. We compare in experiments 1 and 2 the mean thermohaline circulation and oscillatory behavior of a "spherical" ocean sector and a rectangular Cartesian plane, both extending from 20�N to 60�N with similar width at 40�N (5120 km). The forcing is the "linear" surface heat flux, zonally uniform and linearly varying with latitude, but whose area-average is subtracted for the spherical coordinates experiment to insure no drift of the mean basin temperature (the extrema become 40.4 and -49.6 instead of �45 W m-2). In Cartesian coordinates, the Coriolis parameter is set to 2Wsin(latitude) such that the value of b is similar to the one in the spherical case for each latitude.
The convergence of meridians with latitude does not have dramatic consequences on the time-mean circulation, except a slight decrease in the kinetic energy (KE hereafter) and the poleward heat transport compared to the Cartesian experiments. Since the most energetic circulation (western boundary current and its extension in an eastward zonal jet in the polar region) takes place in the polar half of the basin, where the spherical basin is the narrowest, the transports of mass and heat are lowered. There is no common modification in the oscillatory behavior: the periods change by a few years and the index by some 1e-3 �C. The variability pattern (standard deviation of temperature in the upper 250 m) are nearly identical (Fig. 1 and 2).

c. The weak influence of beta
W96 showed the weak influence of the b parameter on the variability and that an f-plane reproduces the oscillatory behavior. We confirm the first result here as a transition between the last Cartesian experiment with b varying realistically and all the following on a constant b-plane centered at 40�N. There are very little differences between the mean states and variability patterns of these two configurations whatever the model used. However, setting b to zero has a noticeable influence, by modifying the extreme values taken by the rotation rate. There is still no agreement in the implied changes in the diagnostics (mean or standard deviation), except a net increase in the mean meridional overturning and consequently, in the mean bottom temperature (because of a reduced residence time at the surface). The change in the standard deviation patterns (Fig. 1 and 3) suggest that the b-effect enhances the variability in the northern regions (where f is the strongest), whereas the variability in the f-plane experiment is spread over all latitudes, with two extrema in the center of the north-west and south-east quarters of the basin.

d. Influence of the Coriolis parameter on Cartesian f-planes
The analysis of the variability in f-plane experiments is simpler because it removes the contribution of b to the vertical velocities, which are then determined through vertical integration of the divergence of the Reynolds stress (in the purely geostrophic model PG0, the vertical velocities thus cancel everywhere except along the boundaries). Then, it is straightforward to compare the influence of the geostrophic and ageostrophic velocities by varying the constant Coriolis parameter (exp. 3 to 7). In all the models, the increase of the rotation rate reduces the part of the ageostrophic terms responsible for the non-zonal circulation. Consequently, KE, meridional overturning and poleward heat transport decrease. Since the surface fluxes have more time to influence the surface temperature, colder waters arise in the northern regions, convect and fill the bottom, while warmer surface temperatures (not mixed down by convection) rise the mean surface temperature. In spite of a less energetic mean circulation, the oscillation amplitude and period increase with f. This suggests that the oscillation relies on the geostrophic balance rather than ageostrophic terms, in agreement with the conclusions regarding the damping influence of increasing the friction (section 3i).

e. The damping effect of the horizontal diffusion
Horizontal tracer diffusion appears to be the most critical damping term for the variability. Whatever the model, forcing and initial state, there is a critical value of KH above which oscillations are damped out (varying between 800 and 2500 m2/s depending on the model used and the horizontal resolution). Estimations of horizontal mixing from buoys trajectories or tracer release experiments are of the order of 1000 m2/s, hence the ocean would be close to a marginally-stable state regarding these oscillations. These are also the values commonly used in models with similar resolution, in which decadal oscillations do not always occur when constant flux forcing is applied. However, higher resolution (marginally eddy-resolving) models, where horizontal mixing is actually carried out by eddies generated by baroclinic instabilities, seems to reproduce the decadal variability as well (Fanning and Weaver 1996). Horizontal diffusion also controls the regularity of the oscillation, its amplitude and period (Fig. 4). In all the models, increasing the horizontal diffusion smoothes the horizontal gradients of temperature and reduces the horizontal velocities (and thus the mean KE). The diffusive part of the poleward heat transport also increases while the deep waters warm and the surface waters cool as the diffusivity increases. This follows since diffusion at high latitudes (where gradients are strong) artificially transports relatively warm waters poleward into the convection regions, so forming relatively warm bottom water. The reduced surface temperatures are due to the reduced poleward heat advection. Regarding the variability, increasing KH lengthens the oscillation period and reduces its amplitude. These results are consistent with the conclusions of HC94. As an alternative to the centered difference advection scheme and explicit horizontal diffusion, we implemented an upstream advection scheme in the horizontal, with no explicit diffusion (the implicit diffusion of such a scheme is U�x/2, reaching 103 m2/s for U=1cm/s and �x=160km). This approach of reducing the explicit uniform diffusion boosts the variability in all the models, in agreement with an average numerical diffusion lower than the previous explicit values. Thus, the oscillation is robust to such a simple numerically-derived velocity-dependent parameterization of horizontal mixing, where the horizontal eddy-diffusivity is increased in regions of strong currents (where more eddies are generated by baroclinic instability). The strong influence of horizontal diffusion is in agreement with an advective mechanism as proposed by Weaver and Sarachik (1991), where density anomalies travel in the polar regions and modify convection, hence the strength of the thermohaline circulation. When diffusion is too strong, these anomalies are diffused before reaching the convection regions, analogous to when too fast a relaxation is applied to the surface thereby damps out any SST deviations. Similarly, Chen and Ghil (1996) show that in their ocean model coupled to an atmospheric Energy Balanced Model, the oscillation amplitude is strongly dependent on the horizontal diffusion coefficient for the atmospheric temperature (playing a role similar to the horizontal eddy-diffusivity in the ocean for the damping of SST anomalies).

f. The driving effect of the vertical diffusion
Vertical diffusion is the engine of the thermohaline circulation, since it is the main source of potential energy (when no wind-induced mixing is considered). Thus, a tight dependency of the strength of the overturning cell is expected. The scaling of the thermocline equations (detailed in Appendix) shows a dependency of the meridional overturning in Kv^(1/2) under fixed-flux boundary conditions (as shown by HC94-we acknowledge a perfect fit to Kv^0.43 for all the models). Thus, under the hypothesis of spatially uniform vertical mixing, the realistic range of variation for the coefficient of diffusivity is narrow in the constraint of realistic overturning. The sensitivity analysis of the oscillation behavior with respect to KV is therefore strongly linked to the sensitivity to the strength of the circulation, and the formulation of the conclusion needs caution: the larger KV, the stronger the circulation and the oscillation, the longer the period. We run all the models for KV varying from 10-5 to 10-3 m2/s. The mean circulation evolves as expected with increasing KV: KE and meridional overturning increase strongly. More vertical mixing of the warm subtropical thermocline induces a cooling of the surface temperature in this area, while the sub-polar temperatures mixed vertically by convection do not change much: thus, the mean surface temperature decreases strongly. The mean bottom temperature, depending principally on the surface water properties in the convection regions, varies very little. Surprisingly, the variability does not vary linearly with KV: all the models show a minimum of the variability for Kv=3e-5 m2/s (PGR0 reached a steady state for KV<5e-5 m2/s). From 3e-5 to 10-3 m2/s, the variability of the temperature field and of all the diagnostics increases faster than linearly, while the oscillation period lengthens (from 17 to 43 years for PGL). These conclusions are in agreement with HC94, although the meridional cell is always much stronger in their study (75 to 180 Sv). In the quest of a mechanism, these experiments are rather disturbing because the doubling of the period between the extreme values of KV does not match the increase by a factor 7 of the meridional overturning, since the advective time scale in the thermohaline loop varies as the inverse of the mass transport: a wave-propagation mechanism, stratification dependent, may be more appropriate.

g. The weak influence of the convective adjustment
Coarse-resolution ocean models cannot resolve the very short time and spatial scales of open-ocean convection. Therefore, it is parameterized independently of the advection-diffusion scheme by instantaneous vertical mixing of unstably stratified grid-boxes at each time-step. It is shown in H96 that this process accounts for 30 to 60% of the potential energy sink (this is roughly the computational cost of this routine as well). Marotzke (1991) shows that the numerical method used for this adjustment influences the stability of the thermohaline circulation under mixed boundary conditions. Cessi (1996) points out a grid-scale instability induced by such instantaneous convective adjustment. In order to assess the robustness of our oscillations to the convection scheme, we ran each model with convective adjustment removed. This is not physically appropriate since static instabilities now remain in the density fields. However, it proves that the variability is not sustained by the convective adjustment.
Without convection, the mean circulation is much more energetic since horizontal density gradients, no longer mixed down, are enhanced (exp. 20 vs. 0, 21 vs. 3): the mean KE, meridional overturning and poleward heat transport are increased. Less cooling reach the bottom waters, which subsequently warms. The surface cannot export its colder waters to depth as efficiently and so cools. PGRW was the only model whose behavior did not change dramatically with or without convection. This is probably related to the high contribution of the KE dissipation in the potential energy depletion (in H96, the role of convection in this model was twice smaller than in the others: downward vertical velocities efficiently drove the coldest water to depth). More relevant to our variability study, the regular oscillations observed in most of the models with convection now become chaotic, although the most energetic time-scales inferred by Fourier analysis remain in the decadal range, usually longer than the period with convection. The amplitude of the oscillation increases strongly for all the diagnostics. These results suggest that convection does not have a driving role in the oscillation, but rather acts as a damping term by mixing negative temperature anomalies downward. However, this is also done by vertical velocities (especially along the boundaries): without convection, these velocities are even more intense and contribute more efficiently to the vertical transport of heat and mass.

h. Influence of the momentum dissipation parameterization and associated boundary conditions
The momentum dissipation and associated boundary conditions significantly influence the mean circulation and deep water characteristics of a thermohaline circulation model forced by restoring boundary conditions (H96). This influence is verified here under fixed flux boundary conditions: the oscillation period (25-34 years) and index (29-56e-3 �C) vary widely depending on the choice of viscosity scheme and the closure for tangential velocities along the boundaries. However, the differences among the time-averaged circulation are not as pronounced under this zonally uniform forcing as under restoring boundary conditions, where the zonally- and meridionally- averaged circulations, as well as deep water properties, were principally related to the upwelling along the western boundary. There, the most intense surface fluxes occurred and had to be compensated for, by upwelling of cold deep water when the dynamics makes it possible. Thus, for poleward heat transport comparable to the restoring runs, the time-averaged meridional overturning varies between 6.8 and 8.7 Sv instead of 7.6-15.6 Sv, but the ordering of the models remains almost the same. Curiously, although the range of variation of the mean bottom temperature is reduced, the ordering is the opposite from the relaxation runs: here, cold bottom temperature is associated with high overturning rate, in opposition to the common sense which would promote colder bottom water due to longer residence time at the surface. The influence of each momentum dissipation scheme or boundary condition on the oscillation features is not straightforward to sort out. Free-slip boundary conditions significantly increase the oscillation amplitude and period, under Rayleigh and Laplacian viscosity (but increase as well the non-linearity of the oscillation, as shown by less sinusoidal time evolution). Such a sensitivity to boundary closures is not well understood yet. i. The damping influence of friction The influence of the friction coefficient is similar, whether the dissipation operator is linear on the horizontal (eH in PGR0) or the vertical (eV in PGS) or harmonic (horizontal AH in PGL, vertical AV in MOM). Increasing the viscosity coefficient always reduces the mean KE. Increasing horizontal viscosity slightly increases the advective poleward heat transport (due to enhanced ageostrophic circulation) and the meridional overturning, while increasing the vertical friction reduces the meridional cell intensity and very weakly the poleward heat transport. The bottom water as well as the surface water warm slightly: this paradox implies a less stratified thermocline, resulting from more vertical mixing induced by enhanced ageostrophic circulation. Increasing the friction coefficient has the opposite effect on the oscillation from increasing the Coriolis parameter in the f-plane experiments: the period shortens and the amplitude decreases (index and all diagnostics standard deviation except the meridional overturning, whose response hesitates between increasing with the mean overturning or decreasing with the variation of the temperature field). We observe here more sensitivity to the mixing of momentum than suggested by HC94 (only the decrease of the mean KE with increased viscosity agrees with their study). The sensitivity of the variability to the vertical mixing of momentum is weak, in spite of the large range of viscosity coefficients tested (AV=0 to 0.07 m2/s in the MOM): the amplitude of KE variations remains lower than 20% while the period changes by less than a year.

j. Influence of the horizontal resolution
The various models do not respond the same way when horizontal resolution is changed. Since some of them (PG0, PGR0 and PGRW) implement "numerical" boundary layers that do not allow a convergence of the circulation when resolution is increased, the mean circulation varies accordingly. Usually, the mean KE, meridional overturning and advective poleward heat transport decrease as well as the mean bottom water temperature with increasing resolution. However, the oscillation behavior is quite robust: the period and the oscillation index vary slightly with finer resolution, except with the linear friction where the variability is reduced for the coarsest grid. The other diagnostics vary widely and unpredictably. In addition to the modifications to the mean circulation depending on the resolution, one must consider the variation of the grid Peclet number (UDx/KH) when resolution, but not the model parameters, is changed: diffusion seems promoted relative to advection when the grid-spacing is reduced. Curiously, HC94 found the oscillation very sensitive to the horizontal resolution in their model: this may be a consequence of the strong variability of the meridional overturning driven by the natural salt boundary condition. Below 100 km grid-spacing, the planetary geostrophic approximation is no longer justified and we observe discrepancies with the primitive equations (PE) dynamics. With higher resolution, PE models seem to induce variability more easily, suggesting that the non-linear terms or time derivatives in the momentum equations help it by providing more inertia to the anomalies. We plan to investigate this issue in a close future.

k. The sensitivity to the vertical discretization
The influence of the vertical discretization is weak when we increase the resolution and the number of levels from the control 15-level run, but quite important if we try to reduce it (exp. 46 to 57). First, this confirms that the choice of vertical discretization used throughout this study is satisfying in terms of the oscillation characteristics presented in Table 2. However, reducing the number of levels and increasing the upper levels thickness leads to strong distortions in the variability and the mean state of the thermohaline circulation. In PGL, the use of 8 levels instead of 15 in the "linear" flux case halves the period and the oscillation index, while in PGR0, it suppresses the variability. When fewer levels are used, the variation of the levels thickness with depth is too important and numerical problems arise when centered differencing is used for the vertical advection: sources of cold water appear at mid-depth along the boundaries where vertical velocities are maximum (Bryan et al. 1975; Weaver and Sarachik 1990; Yin and Fung 1991). In such cases, an upstream scheme was used for the vertical advection: then, the associated implicit vertical diffusion strongly increases the energy of the thermohaline circulation and the overturning rate, as well as the amplitude and period of the oscillations. Thus, even a 2-level discretization reproduces the oscillatory behavior and variability patterns: this proved very useful in the study for the underlying mechanism of the oscillation.

l. Influence of the forcing amplitude
Scaling arguments (detailed in Appendix) show a dependency of the meridional overturning with the amplitude of the surface heat flux to the power 1/4, which agrees with the numerical experiments (58 to 63) where the linear forcing varies from �5 W m-2 to �135 W m-2: The best fit for the overturning as a function of the forcing amplitude is obtained for a regression coefficient of 0.22 in a log-log plot. The mean circulation is increasingly energetic when the heat flux extrema are increased, in KE, poleward heat transport and overturning. The mean surface temperature increases with the stronger warming of the thermocline, while the mean bottom temperature decreases with the more intense cooling at high latitudes. The oscillation amplitude increases regularly with the stronger thermohaline circulation (standard deviations of temperature, KE and overturning), similarly to the effect of increasing the vertical mixing. However, the period of the oscillation evolves oppositely since it becomes shorter when the forcing amplitude increases. The stratification expectedly influences the oscillation period by changing the adjustment waves speed: weaker stratification, associated with higher vertical mixing or lower vertical density gradient (lower forcing amplitude), induces slower adjustment waves that imply longer oscillation periods.

m. Influence of the basin width
The sensitivity of the oscillation to the zonal extent of the basin is an important result in order to give hints for the mechanism of the variability. As pointed out by Frankignoul et al. (1996) in the analysis of the ocean circulation response to stochastic wind forcing, the time required for the adjustment by Rossby waves depends linearly on the basin width. However, by varying the zonal extent of our model basin (reducing the number of grid-points along with the width to keep the resolution constant) from 2,560km to 10,240km, we observe a dependency of the period with the width in only 0.46 power law. As the basin is extended zonally, the time-mean meridional overturning increases (in good agreement with the scaling in Appendix) as well as KE and poleward heat transport, while the mean surface temperature increases to balance the colder deep waters formed further east along the northern boundary. Meanwhile, the oscillation amplitude and period increase as well as the standard deviation of KE and heat transport (it is not systematic for the meridional overturning). Thus, wider basins produce stronger oscillations with longer periods, increasing as the square root of the width instead of proportionally.

4. The role of the boundaries

a. Motivation
W96 and GP96 promote a boundary wave propagation (viscous analogue of Kelvin waves) that would achieve the appropriate time-scale when propagating along the weakly stratified boundaries. However, these two studies do not agree on which boundaries are necessary. While both suggest that the wave does not need to propagate around the whole basin to sustain oscillations, W96 sees the source of the variability in the north east corner of the basin, where the eastward zonal jet extending the western boundary current sinks along the eastern boundary. He presents a convincing description of a perturbation propagating eastward along the northern boundary with the right time-scale. GP96, in a systematic analysis of the role of each boundary, conclude that the western boundary is the only important one for perpetuating the oscillation. In order to shed light on this issue, we analyze the influence of the eastern, northern and southern boundaries successively, and show that none of these is fundamental to the oscillatory behavior. We conclude, in agreement with GP96 but through a totally different analysis and mechanism, that only the western boundary is necessary to generate the variability. However, the intensification of the standard deviation patterns in the eastward zonal jet as an extension of the separated western boundary current suggests that this current is not solely important.

b. The eastern boundary
In order to analyze the influence of the eastern boundary, the basin width is doubled through the inclusion of an eastern 'buffer' zone forced with no surface heat flux (exp. 67). In all models, the consequence of such a change is not dramatic. Figure 5 compares the variability of the surface temperature to the control b-plane experiment with Laplacian friction under the zonally uniform "linear" surface fluxes. The extrema of the standard deviation of surface temperature are slightly shifted, but the oscillation conserves its initial features. In all models, the period lengthens by about 10 years and the amplitude increases compared to the control run, suggesting that the shorter the distance to the eastern boundary, the more damped the oscillation (in agreement with the enhanced variability observed when the basin width is increased in section 3.m). The variability pattern remains almost unchanged in the north-west corner, is the weakest in the south-west corner, and a new area with significant variability appear in the center of the tropical region. The same conclusions apply when more realistic surface heat flux fields are used (diagnosed from the steady-state under restoring boundary conditions for the surface temperature as in exp. 73), although the western intensification of these fluxes traps the variability in the western regions and reduce even more the influence of the eastern boundary (the period increases by 5 to 9 years while the oscillation index does not significantly change).

c. The polar boundary
Following the approach used to examine the role of the eastern boundary, the basin is now extended northward by an equal area with no forcing (exp. 68). The characteristics of the b-plane are conserved, although it leads to unrealistic values of the Coriolis parameter in the northern regions. Nevertheless, no fundamental changes occur in the oscillatory behavior or the mean circulation. The period lengthens by more than 10 years and the oscillation amplitude increases from 30 to 100%. As shown in Fig. 6, the variability pattern is very similar to the standard run, slightly intensified and shifted northward. Obviously, the northern boundary plays a weak role in sustaining the oscillation but influences the period. The meridional overturning adjustment to the north-south pressure gradient is certainly modified depending on the position of the boundary. However, the lack of northern boundary to balance pressure gradients suggest that a viscous boundary wave cannot explain the (almost unchanged) variability in the former boundary region. The stratification of the extended area (N�5�10-4 s-1) should induce time scales much shorter than decadal, would a viscous boundary wave travel in this region. Furthermore, the variability of the temperature field is very weak in this area. Similar conclusions were obtained using the other models and also with more realistic fluxes diagnosed at the end of the restoring runs. The next experiment (#69) was set up with the objective to build a thermocline along both southern and northern boundaries: for this purpose, a symmetric surface forcing was applied to the basin with extended northerly domain, but an f-plane is chosen in order to keep the Coriolis parameter within a geophysical range. A cosine-dependent surface heat flux profile (maximum along the northern and southern boundaries and minimum at mid-latitude) was used. This experiment produced very intense variability with extreme variation of the regions of convection. The oscillation period is longer by 7 years to the one in the same f-plane experiment with a cosine profile on the standard geometry (exp. 70), while the oscillation index doubles. In these cases, we observe a westward propagation of temperature anomalies in the eastward zonal jet. Without a zonal boundary, another mechanism must be found to justify such propagation: recalling classical Rossby waves are not possible since this experiment considers an f-plane, but potential vorticity waves on the mean stratification might be an alternative. d. The western boundary The role of the tropical boundary is investigated in the same manner as the previous eastern and northern ones, by extending the basin of an equal area southward with no surface forcing, the b-plane now spanning from 20�S to 60�N: We do not anticipate a significant dynamical role of this boundary and the numerical simulation (exp. 71) confirms our expectation. Then, only the western boundary remains capable of sustaining the oscillation: GP96, relaxing the temperature field on a very short time scale along each boundary successively, arrived at a similar conclusion. W96 suggested that a periodic basin would not produce oscillations but Marotzke (1990) found decadal variability in the Antarctic circumpolar current, suggesting that boundaries may not be crucial. As a first step towards the understanding of the role of the western boundary, we extended the control basin westward by an equal area with no forcing (exp. 72). The circulation remains strikingly similar to the control run, although the variability is reduced and shifted eastward from the western boundary current, where the cooling fluxes apply. Numerous explanations seems possible at this point, but the intensification of the variability patterns in the north-west quarter (in most of the figures), not along the western boundary current but further eastward in the zonal jet, must be considered. The role of the adjustment of the meridional overturning to the north-south pressure gradient seems important but may just be a consequence of temperature anomalies traveling along the polar boundary. As an alternative to viscous boundary waves and Rossby waves, the adjustment process may imply geostrophic potential vorticity (PV hereafter) waves on the mean stratification, since the northward rising of the thermocline induces a westward propagation of anomalies similar to the b-effect (see section 5d). However, diagnostic calculations of such modes is made difficult in the complex stratification and the associated eastward mean flow. Then, the oscillation could be driven by an instability of the eastward zonal jet occupying the northern half of the basin: Colin de Verdi�re (1986) showed that planetary geostrophic equations are likely to develop baroclinic instability. Another key ingredient, supported by the last experiment, may be the negative surface heat flux that cools the surface as soon as the eastward jet is deflected, since this is the only source of heat in the region. The last feature that might explain the westward propagation of anomalies, even on an f-plane, may be the deep westward current (feeding the deep western boundary current) in the northern regions, as suggested by two-levels experiments.

5. Hints for the mechanism in the transition from steady to oscillatory states

The mechanism invoked by Greatbatch and Zhang (1995) and HC94 involving the zonally-averaged view of an anomaly of density traveling in the mean circulation and influencing the meridional overturning is rather satisfying in the ideal case of a weakly forced and dissipated system, where anomalies circulates endlessly as in the Howard-Malkus loop (Malkus 1972). However, this is not the case on the decadal time-scale (given the large horizontal diffusion used) and moreover, it lacks the fundamental dynamical link between the overturning and the north-south pressure gradient, as well as the consequent time lag between the two, necessary for the adjustment. No variability is observed in a zonally-averaged 2D ocean model forced by constant surface fluxes with traditional "diagnostic" closures for the dynamics (W96 and our own experiments). On the other hand, the experiments of Saravanan and McWilliams (1995) suggest that high frequency noise in the atmospheric forcing is integrated by the ocean in interdecadal time-scales. As pointed out by Cai et al. (1995) in the analysis of the influence of the zonal distribution of the forcing field, the three dimensions are necessary in the adjustment process to build the east-west pressure gradient that controls the zonally-averaged northward transport. This adjustment process that has the same time-scale as the oscillation (spin-up or spin-down of the thermohaline circulation as well) appears as the key mechanism of the variability, as emphasized in recent studies (GP96; W96). We develop only the analysis of the mechanism driving the variability in 3D ocean models, since the 2D ocean response to stochastic forcing might be very different in its roots. a. Shift in the mean state when restoring boundary conditions are changed to fixed fluxes H96 conducted a series of experiments under restoring boundary conditions where the surface temperature was relaxed to a linear profile varying from 25�C at 20�N to 2�C at 60�N, with a relaxation constant of 35 W m-2 K-1 (Haney 1971). We now investigate the transition of the equilibrium states reached upon a switch to the diagnosed surface heat flux boundary condition. The shift to fixed surface heat fluxes does not necessarily lead to variability (only PGRW and PGS do not drift from the steady-state). However in most of the cases, variability occurs spontaneously within a few hundred of years, starting with small amplitudes that grow over several hundred years to a constant amplitude, as the mean basin stratification adjusts to the cooler bottom waters. Curiously, the period of the oscillation does not seem to vary between the weak first oscillations and the final large steady ones (see Fig. 8), and the growth phase of the oscillation amplitude occurs over only a few oscillations. These experiments provide a unique opportunity to see how the steady circulation gets destabilized and evolves into an oscillatory regime. The change in the mean state is almost the same for all the models considered. The mean surface temperature increases along with a decrease in the mean bottom temperature, since cooler minimum temperatures are now achieved at the surface. Potential energy thus decreases, while kinetic energy slightly increases. However, the time averaged available potential energy computed with the method described in H96 usually slightly increases, thus discarding this stability argument to justify the oscillatory state. The time-averaged meridional overturning cell slightly weakens, although the zonal one slightly increases, as does the maximum advective poleward heat transport and the potential energy depletion by convection.

b. Change in the variability patterns as a function of the forcing field
Depending on the model, parameters and forcing, two types of variability appears to take place. The first kind shows a standing temperature anomaly in the north-west region, alternatively positive and negative, as shown in Fig. 9 (a sequence of snapshots of the upper 250 m temperature and velocity anomalies during an oscillation period for PGL under diagnosed heat fluxes). The second kind, that occurs more often on f-planes or under zonally-uniform heat flux, consists of temperature anomalies traveling westward against the eastward zonal-jet in the northern half of the basin. These anomalies induce geostrophically opposite anomalies in their wake by deflecting the mean zonal jet, the only provider of heat in this region of surface cooling. Under diagnosed heat fluxes from the restoring runs, we observe an intensification of the variability in the western regions compared to the variability under zonally-uniform fluxes used in the parameter sensitivity analysis (Fig. 1 and 10d for the b-plane, 3 and 11d for the f-plane, but this conclusion applies for all the other models as well). Since the diagnosed fluxes are also strongly intensified as well in the western regions, especially for the maximum cooling above the warm western boundary current, a relation may exist between the local stability of the circulation and the surface fluxes over the area.

c. The regions with maximum anomaly growth rate
The analysis of the maximum temperature deviations during the first years following the transition to fixed-flux boundary conditions suggests that an instability exists in regions of strong cooling and that the growth rate is maximized where surface fluxes are the lowest. This result is in agreement with the relation suggested previously between variability patterns and surface heat fluxes, although the standard deviation of temperature during the regular oscillations is not necessarily related to the instability mechanism inducing the initial deviation from the steady-state. It also recalls more general conclusions drawn within mixed boundary conditions by Weaver et al. (1991) relating decadal variability to the existence of a minimum in the freshwater fluxes. Similarly, Chen and Ghil (1995, 1996) link the variability to the high latitude regions where strong evaporation or cooling take place, in their sector ocean model forced by mixed boundary conditions or coupled to an atmospheric Energy Balanced Model.
In order to test these ideas further, we looked for the extrema of the surface temperature deviations from the equilibrium state during the first months after the transition to the diagnosed fixed fluxes. We further compared the results from the various models producing different surface fluxes at equilibrium. Most of the extrema for the first months are within one grid-point of the maximum of surface cooling, but starts to be advected by the mean circulation within a few months: Fig. 10 and 11 compare the temperature deviations from the initial equilibrium state to the surface heat flux pattern, diagnosed in the b-plane and f-plane cases. All the models producing variability under diagnosed heat flux show that the seeds of temperature deviation occur in regions of maximum surface cooling. There is a striking correlation between the initial upper-layer temperature deviation and the surface heat flux, although initial temperature deviations are hardly related to the standard deviation of surface temperatures during the regular oscillatory state. Our interpretation is that the instability mechanism shifting the initial state drives the ocean circulation into a resonant mode of the basin, such that the oscillation pattern is the consequence of the long baroclinic waves that achieve the adjustment of the thermohaline circulation. However, since the system is dissipative, an instability mechanism is required during each period to provide the necessary energy for sustaining the regular oscillations.

d. A geostrophic amplification of the western anomalies
Snapshots of the anomalous temperature and velocity fields in the transient phase as well as during the fully developed oscillation (Fig. 9) suggest that the immediate geostrophic adjustment of the horizontal circulation usually acts as an amplifier of temperature anomalies developing in the north-western regions (due to the previous instability mechanism for instance). Since the western boundary current (WBC hereafter) is the provider of heat in the northern half of the basin, cooled at the surface, its variations imply large temperature deviations. Thus, the change in the heat content of the north-west regions is perfectly in phase with the WBC transport, as previously noted by Greatbatch and Zhang (1995). However, since the anomalous temperatures develop a few hundred of kilometers away from the western boundary (where the surface cooling is stronger), the induced geostrophic velocities tend to reinforce the initial WBC deviation. Furthermore, associated with deviations in the western boundary transport, the upwelling along the western wall in the northern regions is modified. This upwelling, sometimes referred to as the "Veronis" effect after Veronis (1975), has a large influence on the temperature in the north-west region by bringing cold water from depth: its role on the whole thermohaline circulation and deep water properties (Boening et al. 1995; H96) suggests that it might also play a role in the variability. Thus, since the upwelling is roughly proportional to the western boundary transport, the induced anomaly along the western boundary is the opposite of the one on the east of the WBC (inducing the transport deviation), such that these dipole increases further the WBC anomaly. Then, different processes interacts for the evolution of temperature anomalies. Geostrophically, because of the mean meridional temperature gradient, a temperature anomaly tends to be advected westward: e.g. a warm anomaly induces southward anomalous motion on its east side, such that cooler waters are advected in the area, while northward motion on its west side brings warmer waters and shifts the anomaly to the west. Such a PV wave propagation is opposed to the eastward mean flow and depending on how large and wide the anomaly grows (controlled by diffusion as well), one effect or the other will lead the advection. For instance, under diagnosed heat fluxes intensified to the west (Fig. 9), the anomalies are quite stationary in the north west corner, while in the control run under zonally uniform heat flux, large anomalies move clockwise across the basin as small opposite anomalies propagate westward along the northern boundary. In addition to the propagation mechanism, given the northern region heat balance between eastward advection of heat and surface cooling, an anomaly deflects the eastward zonal jet in its wake such that an opposite anomaly grows there. When the initial anomaly is trapped along the western boundary, the opposite one, getting large enough, propagates westward and reverses the western boundary current anomaly.

e. The meridional overturning adjustment to the density anomalies
Of course, the geostrophic amplification mechanism previously described cannot hold for long, since it builds its own destruction. On longer time scale, a growing warm anomaly in the north-west corner, increasing the WBC transport, finally reduces the north-south density gradient and modifies the global overturning. However, the adjustment of the meridional cell to the density structure implies geostrophic (long baroclinic Rossby waves and potential vorticity waves) and boundary waves, whose influences are not easily predictable. Therefore, in order to relate the oscillation period to this adjustment process, we further use an indirect method.
The analysis of the response of the thermohaline circulation to small perturbations introduced through the salinity field in the northern regions enlightens the adjustment mechanism. At equilibrium under restoring boundary conditions for surface temperature (salinity uniformly set to 35 psu and not forced), we introduce a 0.05 psu salinity anomaly in the upper 150 m of the region between 50� and 60�N over the longitudinally-centered third of the basin. Salinity is used rather than temperature because it does not directly change the surface fluxes under restoring boundary conditions and allows us to differentiate the introduced perturbation from the dynamically induced one. A linear equation of state relates the density anomaly to the temperature (a=2e-4 K-1) and to the salinity (b=8e-4 psu-1). Then, we compare the adjustment of the equilibrium state to the weak salinity anomaly under the restoring boundary conditions to the associated diagnosed surface flux boundary conditions. These experiments show a similarity in the response during the first years, but a different evolution of the perturbations afterward (Fig. 12): The anomalous temperatures are damped out within a few cycles under restoring boundary conditions, while they amplify under fixed fluxes and reach the full amplitude of the oscillatory state within a few hundred years, the time for the mean basin temperature to adjust as well. Except in PGRW, we observe an overshoot of the meridional overturning in response to the perturbation under both boundary conditions. But these experiments suggest mainly that the period of the oscillation is set by the mean stratification and circulation, independently of the perturbation amplitude: even the infinitesimal first cycles on Fig. 12 show the same period as the fully developped oscillations, such that the temperature anomalies amplitude do not influence their dynamics. This is noteworthy since most of the regular oscillations induce large changes in the upper temperature and circulation, such that non-linear interactions might be expected as compared to the small perturbation experiment done here. Temperature anomalies associated with geostrophic circulation anomalies are seen traveling westward against the eastward zonal jet in the northern part of the basin and crossing the basin within a few years before reaching the western boundary current: this looks more like geostrophic PV waves than viscous boundary waves, since there is no intensification close to the coast.
This behavior is observed in f-plane simulations as well, thus discarding the classical baroclinic Rossby wave explanation. Nevertheless, potential vorticity waves on the mean circulation-stratification might provide an alternative mechanism for the westward propagation, the oscillation time scale and for shaping up the oscillation pattern. The two major features one must consider in this argument are the poleward rise of the thermocline depth (in analogy with a topographic beta-effect) and the associated geostrophic eastward jet in the upper layers.

f. A simple box-model analogy for the oscillation
Drastically reducing the degrees of freedom of ocean models have proved successful in the past, for the understanding of the multiple equilibria of the thermohaline circulation for instance (Stommel 1961). Box models are the most simple architecture to investigate process studies and we try in the following to brush a picture of the decadal oscillation in such an idealized formalism. The analysis of these oscillations in terms of zonal averages is complicated by the phase-shift of temperature deviations along the same latitude-circle, such that zonally-averaged temperatures varies much less than expected from local deviations. Thus, the definition of mean temperature for the traditional polar (T1) and tropical (T2) boxes in these models filters most of the variability. However, it allows us to look at the correlation between meridional overturning (O) and north-south density gradients (proportional to DT=T2-T1). We define a polar and a tropical box of equal size, separated by the 40�N parallel and 400 m deep. In Fig. 13, we plot the time evolution of the mean temperatures in these boxes along with their difference, the meridional overturning, the WBC transport in the upper 400 m and their time derivatives, for the low-resolution run under linear heat flux (exp. 42). In contradiction with what is usually stated in box and zonally-averaged ocean models, meridional overturning anomalies and north-south density gradient anomalies are closer to quadrature (correlation coefficient~0.95) than in-phase (c~0.03). Other models give similar results although the time evolution of the meridional overturning is less sinusoidal at higher resolution (see Fig. 8 vs. 13). This confirms that on decadal time-scales, the adjustment process cannot be neglected and the equilibrium overturning is never reached. Consequently, a simple box-model analogy representing the decadal oscillations is formulated as follows, with a parameterization of the adjustment process such that the time-derivative of the meridional overturning is linearly related to the north-south density gradient anomaly. We consider a thermocline box (volume V�2000km�4000km�400m, temperature T=12�C), warmed by a constant heat flux (A*QH, A=2000km�4000km, QH=35 W m-2), and a cold box (volume V0>>V, temperature T0=5�C) including the polar and deep waters, cooled by the opposite heat flux, such that the mean basin temperature remains constant (Fig. 14). The mass exchange (O�10 Sv) between these boxes represents the traditional overturning cell, moving thermocline waters poleward through the western boundary current, balanced by upwelling across the base of the thermocline (in the interior and/or along the western boundary). Thus, the equations controlling the temperature of the boxes are:
Decomposing T, T0 and O in a time-average and deviation, respectively noted by an overbar and a prime, the time-averaged equations become: while the first order linearization for the deviation from the time-average gives: (E1) recalling that the mean basin temperature is constant:. From the previous paragraph, the rate of change of the overturning deviation is chosen proportional to the south-north temperature anomaly gradient to parameterize crudely the time-lag between anomalous pressure anomalies and overturning adjustment (the proportionality constant k is found around 1 Sv year-1 K-1 in most of the numerical experiments): (E2) In a first approximation mostly verified in the experiments forced by diagnosed heat flux (Fig. 8), we suppose that (this is barely achieved in fig. 13 by a factor of 3), such that the second term in the right hand side of (E1) can be neglected. Then, the time evolution of T' is simply sinusoidal with period �23 years, with =7�C. However, when solving the full system (E1, E2), the previously neglected term in (E1) leads to an e-folding time-scale similar to the period: �19 years. Obviously, we lack a driving term representing the instability processes that supply energy to the perturbation during each period to sustain the oscillation. Nevertheless, this simple model illustrates how crucial for the oscillation is the time lag between the overturning and the north-south density gradient, through the westward propagation of various waves (geostrophic, Rossby, PV or viscous boundary) across the basin. If O is taken in phase with (T'-T0'), no periodic solution can be found. Finally, the parameterization makes sense if one imagines the thermohaline cell forced by differential gravity forces, like in the Welander's circulation tube of Malkus (1972) where the rate of change of the transport is related to the torque of the weight applied to the non-uniform fluid density.

6. Robustness of the variability to various parameterizations and forcing

The range of forcing leading to variability in the literature suggests that at least one variable of the density equation requires more freedom than allowed by the usual short time scales used in restoring boundary conditions. Applying the idea of Schopf (1983), we reduced the implicit infinite heat capacity of the atmosphere under restoring boundary conditions by reducing the relaxation constant (lengthening the restoring time scale). At the same time the specified atmospheric temperature amplitude was increased such that the thermohaline circulation kept a reasonable intensity. Greatbatch and Zhang (1995) carried out similar experiments and reproduced decadal oscillations. With the same parameters as the control run described in section 2, we observed regular decadal oscillations for thousands of years, with amplitude and periods similar to the control run, for restoring temperatures varying from 29�C at 20�N to -2�C at 60�N and a relaxation constant of 10 W m-2 K-1, as suggested by Seager et al. (1995). The standard deviation of the upper 250 m temperature is shown in Fig. 15 (exp. 79) and looks very similar to previous variability patterns. The time-averaged surface fluxes vary between 50 W m-2 in the north-west corner to -30 W m-2 along the tropical boundary. Other cases where we observed variability under restoring boundary conditions with traditional relaxation constant (35 W m-2 K-1) involved those with increased vertical mixing. When the intensity of the meridional overturning becomes strong enough, decadal oscillations occur spontaneously: although in this case the damping of anomalies did not vary, it is the amplitude of the overturning (non-linearity of the western boundary current or eastward zonal jet) that triggers the oscillations.
The potential role of the upwelling along the western boundary in sustaining the oscillation is invoked in the description of the geostrophic processes amplifying temperature anomalies. This is also supported by the reluctance of PGRW (where vertical velocities are reduced along the boundaries by allowing horizontal recirculations in tangential velocities along the coasts) to sustain variability under fixed fluxes diagnosed at the equilibrium under restoring boundary conditions (also if horizontal diffusion is increased to 1000 m2/s under "linear" heat fluxes, PGRW is the only one to remain in a steady state). Therefore, in order to investigate the sensitivity of the oscillatory behavior to diapycnal mixing (through cross-isopycnal vertical velocities along the boundaries), we compared various experiments with the GFDL Modular Ocean Model. The geometry, parameters and forcing remained the same as in the control run, except that now the mixing tensor was rotated along isopycnal coordinates (the background vertical diffusivity was not changed from the control run). The oscillation proved robust to this parameterization: the periods becomes irregular (between 16 and 25 years) and the amplitude increases. We then included the parameterization for mesoscale eddy-induced tracer transport (Gent and McWilliams 1990). Although the amplitude of the oscillation was reduced and the period shortened (16 years instead of 25), the variability was sustained with patterns similar to the control run.
The addition of a steady wind stress, forcing a barotropic streamfunction that interacts non-linearly with the thermohaline circulation, was also carried out. Various profiles of zonally-uniform wind stress were tested, inducing one or two gyres in the standard flat-bottomed basin. In all the cases, the oscillation survives, with shift in the variability patterns, slight changes in the period and increased amplitude. The wind forcing, as implemented here, seems not to influence the mechanism driving decadal variability.
Finally, we introduced seasonal variations in the forcing field (similar to North Atlantic annual cycle for the zonal mean net downward heat flux) to test the robustness of the decadal variability to continuously varying forcing and circulation. This had no significant effect on the oscillation, although the linear equation of state, the lack of wind-stress and freshwater forcing (and their seasonal variations) might contribute to these limited changes. The seasonal cycle, whose amplitude remains smaller than the decadal one, seems not to interact at all with the longer scales of decadal oscillations.

7. Conclusion

We have analyzed in this study decadal oscillations forced by fixed surface fluxes (or restoring boundary conditions with realistically long restoring time scales) in box-geometry ocean basins, with planetary geostrophic and primitive equations models at coarse resolution. Through a systematic comparison of the results from several models with different parameterizations of momentum dissipation and associated boundary conditions and through an extensive parameter sensitivity analysis, we were able to determine which processes are critical to the variability.
Spherical or Cartesian coordinates, as well as the variation of b with latitude or a constant mid-latitude value, were shown to have no influence on the oscillation. However, f-planes experiments suggest a leading role of the geostrophic balance, since increasing the Coriolis parameter increases the oscillation amplitude and period although it weakens the mean circulation. Vertical diffusion has a driving role on the variability, while horizontal diffusion has a crucial damping role. Within the range of uncertainty in the actual ocean mixing, the variability is sustained in all the models or in none. A numerical scheme relating the horizontal mixing to local velocities, as well as isopycnal mixing with or without the eddy-induced tracer transport parameterization of Gent and McWilliams (1990), does not profoundly influence the oscillatory behavior. Convective adjustment, well-known to produce grid-scale noise and influence the stability of the thermohaline circulation, is not necessary to sustain the oscillations (but to the removal of static instabilities).
Momentum dissipation parameterizations and associated boundary conditions do have a significant influence on the oscillation period and amplitude-surprisingly, since they modify the velocities along the boundaries but not the geostrophic balance in the interior. However, the variability is reproduced in all these models. The friction coefficients have a damping effect on the oscillations, in agreement with the driving role of the Coriolis parameter. Thus, decadal oscillations stand out as a robust geostrophic feature of thermohaline circulation models forced by constant surface fluxes in idealized geometry: horizontal and vertical resolutions are not crucial but may modify the variability if not fine enough-the variability occurs more easily at higher resolution.
We further investigated the role of the boundaries in the propagation of viscous boundary waves, as proposed in previous studies. Our results suggest that the propagation of boundary waves along weakly stratified boundary is not always appropriate to explain the variability and time-scale of the oscillation. The intensification of the standard deviation of the upper layers temperature in the north-west quarter of the basin, in the eastward zonal jet continuing the western boundary current where surface cooling is the strongest, suggests instability processes in this area.
We finally look for hints for the mechanism in the analysis of the transition between steady-state (under restoring boundary conditions) and oscillatory regime under the fixed fluxes diagnosed at the equilibrium of the restoring run. The time-averaged state does vary between the two regimes, since cooler deep water is formed during the oscillations (however, available potential energy is not systematically lower in the final state). Initial deviation from the steady-state are shown to be related to the regions of strongest surface cooling, where a purely geostrophic amplification mechanism of temperature anomalies in the vicinity of the western boundary appears to take place. The oscillation period seems related to the meridional overturning adjustment to density anomalies. Finally, a simple box-model is proposed that captures the phase-shift between the meridional overturning anomalies and anomalous north-south temperature gradients and reproduces the oscillatory behavior.

Acknowledgments. This research was funded through a scholarship awarded to TH from the French Ministry of Education and strategic and operating grants awarded to AJW from NSERC, CICS and the NOAA Scripps-Lamont Consortium on the Ocean's Role in Climate.


APPENDIX: Geostrophic scaling of overturning streamfunction under constant flux boundary condition

Huang and Chou (1994) developed a geostrophic scaling for the ocean circulation driven by constant freshwater flux, applied by Weaver and Garrett (1994) to the buoyancy flux. We briefly recall this scaling to allow a comparison with the numerical experiments. Let U and V be characteristic horizontal geostrophic scales, W a characteristic vertical velocity, D a characteristic thermocline depth scale, LX and LY the zonal and latitudinal scales of the ocean basin, f the Coriolis parameter, DrX and DrY characteristic density gradients in each direction (a priori different), and QB the characteristic surface buoyancy flux. From hydrostatic and geostrophic balances, the thermal wind equation gives: or .
The continuity equation relates U, V and W through the typical length scales: . Given a constant vertical diffusivity KV, the Munk balance for the vertical advective-diffusive tracer balance across the thermocline imposes: KV = W D. Then , the density gradient can be related to the surface fluxes via the boundary condition: , i.e. or via the horizontal tracer advection in the thermocline: , which shows that DrX=DrY=Dr. Thus, the thermocline depth and the meridional overturning can be deduced from the parameters: , .
The numerical experiments with various models, considering the variations of each parameter X, give the following power laws, as the regression coefficient of Y function of X in log-log plot:

POWER LAW Parameter [range] Scaling PGL PG0 PGR0 PGRW mean f [2Wsin(20�)-2Wsin(60�)] -1/4 -0.56 -0.35 -0.57 -1.07 -0.64 KV [10-5-10-3 m2/s] 1/2 0.43 0.49 0.43 0.45 0.45 QB [5-135 W/m2] 1/4 0.22 0.32 0.17 0.22 0.23 LX [2,560-10,240 km] 3/4 0.69 0.60 0.73 0.85 0.72

Most of the numerical values match reasonably the scaling except for the Coriolis parameter, a rather disturbing result since the geostrophic equations are precisely the basis for the scaling. The different dynamics between beta- and f-plane might explain the discrepancy.
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, K., S. Manabe, and R. C. Pacanowski, 1975: A global ocean-atmosphere climate model. Part II: The oceanic circulation. J. Phys. Oceanogr., 5, 30-46. Cai, W., R. J. Greatbatch, and S. Zhang, 1995: Interdecadal variability in an ocean model driven by a small, zonal redistribution of the surface buoyancy flux. J. Phys. Oceanogr., 25, 1998-2010. Cessi, P., 1996: Grid-scale instability of convective-adjustment schemes. J. Mar. Res., 54, 407-420. Chen, F., and M. Ghil, 1995: Interdecadal variability of the thermohaline circulation and high-latitude surface fluxes. J. Phys. Oceanogr., 25, 2547-2568. Chen, F., and M. Ghil, 1996: Interdecadal variability in a hybrid coupled ocean-atmosphere model. J. Phys. Oceanogr., 26, 1561-1578. 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. Delworth, T., S. Manabe, and R. J. Stouffer, 1993: Interdecadal variations of the thermohaline circulation in a coupled ocean-atmosphere model. J. Climate, 6, 1993-2011. Deser, C., and M. L. Blackmon, 1993: Surface climate variations over the North Atlantic Ocean during winter: 1900-1989. J. Climate, 6, 1743-1753. Fanning, A. F., and A. J. Weaver, 1996: A horizontal resolution and parameter sensitivity study of heat transport in an idealized coupled climate model. J. Climate, submitted. Frankignoul, C., P. Muller, and E. Zorita, 1996: A simple model of the decadal response of the ocean to stochastic wind forcing. J. Phys. Oceanogr., in press. Gent, P. R., and J. C. McWilliams, 1990: Isopycnal mixing in ocean circulation models. J. Phys. Oceanogr., 20, 150-155. Greatbatch, R. J., and K. A. Peterson, 1996: Interdecadal variability and oceanic thermohaline adjustment. J. Geophys. Res. (Oceans), 101, 20467-20482. Greatbatch, R. J., and S. Zhang, 1995: An interdecadal oscillation in an idealized ocean basin forced by constant heat flux. J. Climate, 8, 81-91. Griffies, S. M., and K. Bryan, 1997: Predictability of North Atlantic multidecadal climate variability. Science, 275, 181-184. Haney, R., 1971: Surface thermal boundary condition for ocean circulation models J. Phys. Oceanogr., 1, 241-248. Hansen, D. V., and H. F. Bezdek, 1996: On the nature of decadal anomalies in North Atlantic sea surface temperature. J. Geophys. Res., 101, 8749-8758. Hasselmann, K., 1982: An ocean model for climate variability studies. Prog. Oceanog., 11, 69-92. Huang, R. X., and R. L. Chou, 1994: Parameter sensitivity of the saline circulation. Climate Dyn., 9, 391-409. Huang, R. X., 1993: Real freshwater flux as a natural boundary condition for the salinity balance and thermohaline circulation forced by evaporation and precipitation. J. Phys. Oceanogr., 23, 2428-2446. Huck, T., A. J. Weaver, and A. Colin de Verdi�re, 1996: The effect of different parameterizations and boundary conditions applied to the momentum equation in coarse-resolution thermohaline circulation models. J. Phys. Oceanogr., submitted. Kushnir, Y., 1994: Interdecadal variations in North Atlantic Sea Surface Temperature and associated atmospheric conditions. J. Climate, 7, 141-157. Latif, M., and T. P. Barnett, 1996: Decadal variability over the North Pacific and North America: Dynamics and predictability. J. Climate, in press. Lenderink, G., and R. J. Haarsma, 1996: On the mechanism of decadal oscillations in a coarse resolution ocean model. J. Phys. Oceanogr., submitted. Malkus, W. V. R., 1972: Non-periodic convection at high and low Prandtl number. M�moires Soci�t� Royale des Sciences de Li�ge, 6, IV, 125-128. Marotzke, J., 1990: Instabilities and multiple equilibria of the thermohaline circulation. Ph.D. thesis dissertation, Institut fur Meereskunde, Kiel, 126 pp. Marotzke, J., 1991: Influence of convective adjustment on the stability of the thermohaline circulation. J. Phys. Oceanogr., 21, 903-907. Pacanowski, R. C., 1995: MOM 2 Documentation, User's Guide and Reference Manual. GFDL Ocean Technical Report #3. Phillips, N. A., 1963: Geostrophic motion. Rev. Geophys. Space Phys., 1, 123-176. Power, S. B., and R. Kleeman, 1994: Surface heat flux parameterizations and the response of ocean general circulation models to high-latitude freshening. Tellus, 46A, 86-95. Rahmstorf, S., 1993: A fast and complete convection scheme for ocean models. Ocean Modelling, 101 (unpublished manuscript). Rahmstorf, S., and J. Willebrand, 1995: The role of temperature feedback in stabilizing the thermohaline circulation. J. Phys. Oceanogr., 25, 787-805. Reverdin, G., D. Cayan, and Y. Kushnir, 1997: Decadal variability of hydrography in the upper northern North Atlantic 1948-1990. J. Geophys. Res., 102, 8505-8531. 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. Saravanan, R., and J. C. McWilliams, 1995: Multiple equilibria, natural variability, and climate transitions in an idealized ocean-atmosphere model. J. Climate, 8, 2296-2323. Schopf, P. S., 1983: On equatorial waves and El Nino. II: Effects of air-sea thermal coupling. J. Phys. Oceanogr., 13, 1878-1893. Seager, R., Y. Kushnir, and M. A. Cane, 1995: On heat flux boundary conditions for ocean models. J. Phys. Oceanogr., 25, 3219-3230. Stommel, H., 1961: Thermohaline convection with two stable regimes of flow. Tellus, 13, 224-230. Timmermann, A., M. Latif, R. Voss, and A. Grotzner, 1996: North Atlantic interdecadal variability: A coupled air-sea mode. J. Climate, submitted. 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 C. J. R. Garret, 1994: On the sensitivity of the thermohaline circulation in OGCM to the vertical eddy diffusivity. unpublished manuscript. Weaver, A. J., and E. S. Sarachik, 1990: On the importance of vertical resolution in certain OGCM. J. Phys. Oceanogr., 20, 600-609. 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. Weaver, A. J., E. S. Sarachik, and J. Marotzke, 1991: Freshwater flux forcing of decadal and interdecadal oceanic variability. Nature, 353, 836-838. Winton, M., 1993: Numerical Investigations of Steady and Oscillating Thermohaline Circulations. Ph.D. thesis dissertation, University of Washington, 155 pp. 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., 1996: The role of horizontal boundaries in parameter sensitivity and decadal-scale variability of coarse-resolution ocean general circulation models. J. Phys. Oceanogr., 26, 2 289-304. Wohlleben, T., and A. J. Weaver, 1995: Interdecadal climate variability in the subpolar North Atlantic. Clim. Dyn., 11, 459-467. Yin, F. L., and I. Fung, 1991: Net diffusivity in ocean general circulation models with nonuniform grids. J. Geophys. Res., 96, 10, 773-776. Yin, F. L., and E. Sarachik, 1995: Interdecadal thermohaline oscillations in a sector ocean general circulation model: advective and convective processes. J. Phys. Oceanogr., 25, 2465-2484. Zhang, S., C. A. Lin, and R. J. Greatbatch, 1992: A thermocline model for ocean-climate studies. J. Mar. Res., 50, 99-124. Zhang, S., R. J. Greatbatch, and C. A. Lin, 1993: A reexamination of the polar halocline catastrophe and implications for coupled ocean-atmosphere modeling. J. Phys. Oceanogr., 23, 287-299. TABLE CAPTIONS Table 1. Experiments summary. Geometry: f states for f-plane while b states for b-plane, followed by centered latitude. Longitude extent: C=5,120km, W=10,240km, N=2,560km. Latitude extent: C=20-60�N, W=20-100�N, X=-20:100�N. Number of gridpoints: horizontally, C=32*28. Forcing: L for zonally-uniform fixed surface fluxes linearly varying with latitude (Cos for cosine profile) followed by extrema (W/m2), R for restoring boundary conditions for surface temperature followed by extrema (�C)-the restoring constant is given in the specifications-, D for diagnosed fixed fluxes at the equilibrium of the associated restoring experiment. KH, KV: horizontal and vertical tracer diffusivities. AH: horizontal viscosity, the vertical viscosity is zero except in the MOM runs when it is stated in the specifications. Table 2. Variability diagnostics for the experiments detailed in Table 1. Oscillation period (year, a ~ prefix means the oscillation is not regular), oscillation index (10-3 �C, as defined in section 3), kinetic energy density (time-average and standard deviation, in 0.1 kg m-1 s-2), meridional overturning (time-average and standard deviation, in Sv=106 m3/s), maximum advective poleward heat transport (time-average and standard deviation, in PW), area-averaged surface (0-50m) and bottom (3950-4500m) temperatures (time-average and standard deviations, in �C). Note that in the fixed flux experiments only, the mean basin temperature is set to 4�C by initial conditions. Figure captions FIG. 1. Standard deviation of the upper 250 m temperature (�C) for the control run (exp. 0): PGL under the "linear" surface heat flux (zonally uniform and linearly varying with latitude from 45 W/m2 at 20�N to -45 W/m2 at 60�N) on a Cartesian coordinates b-plane (centered at 40�N). Regions of variability greater than the area mean plus standard deviation are shaded. The time-averaged velocity field is superimposed (the scale is 3 cm/s per degree). FIG. 2. Standard deviation of the upper 250 m temperature (�C) for PGL under the "linear" surface heat flux in spherical coordinates (exp. 1). FIG. 3. Standard deviation of the upper 250 m temperature (�C) for PGL under the "linear" surface heat flux on a Cartesian coordinates f-plane centered at 40�N (exp. 4). FIG. 4. The influence of horizontal diffusivity on the oscillatory behavior in the standard PGL experiment: Kinetic energy density as a function of time for KH=350 (exp. 8, dashed), 700 (control run, solid), 1500 (exp. 11, dash-dotted) and 3000 (exp. 14, solid and constant) m2/s. FIG. 5. Same as Fig. 1 but for a basin extended eastward by an equal area with no forcing (67). FIG. 6. Same as Fig. 1 but for a basin extended northward by an equal area with no forcing (68). FIG. 7. Same as Fig. 1 but for an f-plane basin extended northward by an equal area, forced by surface heat flux with a cosine profile symmetric relatively to the previous northern boundary (69). FIG. 8. Transition from the steady-state under restoring boundary conditions to the regular oscillations under fixed fluxes diagnosed at the equilibrium of the restoring run (Fig. 10a): meridional overturning, area-averaged surface (0-50m) and bottom (3950-4500m) temperatures, maximum advective poleward heat transport, potential energy depletion by convection and kinetic energy density as functions of time (after the shift from restoring to fixed flux boundary condition). FIG. 9. Snapshots of the anomalies of temperature and velocity in the upper 250 m during an oscillation of PGL (every 2.3 years), forced by the heat fluxes diagnosed at the equilibrium under restoring boundary conditions (Fig. 10a, exp. 73). The scale for velocities is 1 cm/s per degree. Cold anomalies (<-2�C) are in black while warm anomalies (>2�C) are in white. FIG. 10. Transition from restoring boundary conditions to fixed fluxes for PGL on the b-plane: (a) Surface heat fluxes (W/m2) diagnosed at the equilibrium under restoring boundary conditions (exp. 73); (b) Standard deviation of the upper 250 m temperature (�C) in the regular oscillatory regime; (c) and (d) Upper 250 m temperature (10-3 �C) and velocity (10-4 m/s per gridspacing) deviations, one and three month after the shift from restoring to diagnosed fixed flux boundary condition. The maximun cooling and deviations regions are shown by a '+'. FIG. 11. Same as Fig. 10 but on the f-plane (exp. 77 and 78). FIG. 12. Meridional overturning (Sv) as a function of time under restoring boundary conditions (dash-dotted) vs. fixed heat fluxes diagnosed from the equilibrium of the restoring run (dashed), after the introduction of a 0.05 psu salinity anomaly in the centered north of the basin (exp. 75). The responses to the opposite anomaly are plotted as well (dotted). The initial state is the equilibrium state under restoring boundary conditions. The regular oscillations regime under diagnosed heat fluxes has a period of 23 years, very similar to the initial cycles period. Note the damping effect of the restoring boundary conditions on temperature anomalies, compared to the amplification under fixed fluxes. FIG. 13. Mean temperature in the upper 400 m for the north (T1) and south (T2) halves of the basin (�C) and difference, western boundary current transport in the upper 400 m (GS, dash-dotted) and maximum of the meridional overturning streamfunction (O, in Sv), d(T2-T1)/dt (in K/year), dO/dt (in Sv/year) as functions of time during a few oscillations in the low resolution PGL run forced by the "linear" heat flux (exp. 42). The linear regression coefficient is 1.38 (1.0) Sv year-1 K-1 between dO/dt (dGS/dt) and T2-T1, with a correlation coefficient of 0.92 (0.80), while it is -0.04 K year-1 Sv-1 between d(T2-T1)/dt and O or GS, with correlation coefficients of respectively -0.97 and -0.79. The slow growth but fast breakdown of the meridional overturning seen at higher resolution (fig. 8) is avoided here, thanks to the low resolution. FIG. 14. A simple box-model for the oscillation. FIG. 15. Same as Fig. 1 for PGL, under restoring boundary conditions for the surface temperatures to a linear profile varying from 29�C at 20�N to -2�C at 60�N, with a 10 W m-2 K-1 relaxation constant . Table 1. Experiments summary. exp # Geometry Lon ext Lat ext NX�NY NZ Forcing [W m-2] KH [m2/s] KV [m2/s] AH [m2/s] Specifications 0 b 40�N C C C 15 L �45 700 1e-4 1.5e5 Control run 1 Spherical 40� C C 15 L 40/-50 700 1e-4 1.5e5 2 Cartesian C C C 15 L �45 700 1e-4 1.5e5 3 f 40�N C C C 15 L �45 700 1e-4 1.5e5 4 f 20�N C C C 15 L �45 700 1e-4 1.5e5 5 f 30�N C C C 15 L �45 700 1e-4 1.5e5 6 f 50�N C C C 15 L �45 700 1e-4 1.5e5 7 f 60�N C C C 15 L �45 700 1e-4 1.5e5 8 b 40�N C C C 15 L �45 350 1e-4 1.5e5 9 b 40�N C C C 15 L �45 500 1e-4 1.5e5 10 b 40�N C C C 15 L �45 1000 1e-4 1.5e5 11 b 40�N C C C 15 L �45 1500 1e-4 1.5e5 12 b 40�N C C C 15 L �45 2000 1e-4 1.5e5 13 b 40�N C C C 15 L �45 2500 1e-4 1.5e5 14 b 40�N C C C 15 L �45 3000 1e-4 1.5e5 15 b 40�N C C C 15 L �45 0 1e-4 1.5e5 Upstream horizontal advection scheme 16 b 40�N C C C 15 L �45 700 1e-5 1.5e5 17 b 40�N C C C 15 L �45 700 3e-5 1.5e5 18 b 40�N C C C 15 L �45 700 3e-4 1.5e5 19 b 40�N C C C 15 L �45 700 1e-3 1.5e5 20 b 40�N C C C 15 L �45 700 1e-4 1.5e5 No convective adjustment 21 f 40�N C C C 15 L �45 700 1e-4 1.5e5 No convective adjustment 22 b 40�N C C C 15 L �45 700 1e-4 1.5e5 MOM 23 b 40�N C C C 15 L �45 700 1e-4 1.5e5 PGL + free-slip b.c. 24 b 40�N C C C 15 L �45 700 1e-4 0.0 PG0 (purely-geostrophic) 25 b 40�N C C C 15 L �45 700 1e-4 0.0 PG0 + free-slip b.c. 26 b 40�N C C C 15 L �45 700 1e-4 0.0 PGR0 (linear friction) 27 b 40�N C C C 15 L �45 700 1e-4 0.0 PGR + free-slip b.c. 28 b 40�N C C C 15 L �45 700 1e-4 0.0 PGRW no-normal-flow b.c. 29 b 40�N C C C 15 L �45 700 1e-4 0.0 PGS (Salmon,86,90) eV=0.03 30 b 40�N C C C 15 L �45 700 1e-4 0.0 PGS (Salmon,86,90) eV=0.3 31 b 40�N C C C 15 L �45 700 1e-4 1.0e3 32 b 40�N C C C 15 L �45 700 1e-4 1.0e4 33 b 40�N C C C 15 L �45 700 1e-4 5.0e4 34 b 40�N C C C 15 L �45 700 1e-4 1.0e5 35 b 40�N C C C 15 L �45 700 1e-4 5.0e5 36 b 40�N C C C 15 L �45 700 1e-4 1.0e6 37 b 40�N C C C 15 L �45 700 1e-4 1.5e5 MOM AV=1e-4 m2/s 38 b 40�N C C C 15 L �45 700 1e-4 1.5e5 MOM AV=1e-3 m2/s 39 b 40�N C C C 15 L �45 700 1e-4 1.5e5 MOM AV=1e-2 m2/s 40 b 40�N C C C 15 L �45 700 1e-4 1.5e5 MOM AV=5e-2 m2/s 41 b 40�N C C C 15 L �45 700 1e-4 1.5e5 MOM AV=7e-2 m2/s 42 b 40�N C C 16�14 15 L �45 700 1e-4 1.5e5 Low resolution 43 b 40�N C C 24�21 15 L �45 700 1e-4 1.5e5 Intermediate low resolution 44 b 40�N C C 48�42 15 L �45 700 1e-4 1.5e5 Intermediate high resolution 45 b 40�N C C 64�56 15 L �45 700 1e-4 1.5e5 High resolution 46 b 40�N C C C 90 L �45 700 1e-4 1.5e5 dz=50m�90 47 b 40�N C C C 30 L �45 700 1e-4 1.5e5 dz=25x6,(50:25:250) �2,275�6 48 b 40�N C C C 20 L �45 700 1e-4 1.5e5 dz=25�4,50�4,100:50:450,500�4 49 b 40�N C C C 18 L �45 700 1e-4 1.5e5 dz=50�6,100:50:450,500�4 50 b 40�N C C C 15 L �45 700 1e-4 1.5e5 dz=100�4,150:50:450,500�4 51 b 40�N C C C 8 L �45 700 1e-4 1.5e5 dz=50,100,250,450,650,850,1100�2 52 b 40�N C C C 8 L �45 700 1e-4 1.5e5 dz=100�2,200:200:1000,1300 53 b 40�N C C C 8 L �45 700 1e-4 1.5e5 dz=100�2,200:200:1000,1300 Ups.vert.adv. 54 b 40�N C C C 4 L �45 700 1e-4 1.5e5 dz=100,300,1000,3100 Upstream vert. adv. 55 b 40�N C C C 4 L �45 700 1e-4 1.5e5 dz=200,400,800,3100 Upstream vert. adv. 56 b 40�N C C C 2 L �45 700 1e-4 1.5e5 dz=250,4250 Upstream vert. adv. 57 b 40�N C C C 2 L �45 700 1e-4 1.5e5 dz=500,4000 Upstream vert. adv. 58 b 40�N C C C 15 L �5 700 1e-4 1.5e5 59 b 40�N C C C 15 L �15 700 1e-4 1.5e5 60 b 40�N C C C 15 L �30 700 1e-4 1.5e5 61 b 40�N C C C 15 L �60 700 1e-4 1.5e5 62 b 40�N C C C 15 L �90 700 1e-4 1.5e5 63 b 40�N C C C 15 L �135 700 1e-4 1.5e5 64 b 40�N N C 16�28 15 L �45 700 1e-4 1.5e5 basin 50% narrower 65 b 40�N W C 48�28 15 L �45 700 1e-4 1.5e5 basin 50% wider 66 b 40�N W C 64�28 15 L �45 700 1e-4 1.5e5 basin twice wider (Pacific) 67 b 40�N W C 64�28 15 L �45 700 1e-4 1.5e5 Zero heat flux on eastern half 68 b 40�N C W 32�56 15 L �45 700 1e-4 1.5e5 Zero heat flux on northern half 69 f 40�N C W 32�56 15 Cos �45 700 1e-4 1.5e5 Forcing=cosine symmetric / 40�N 70 f 40�N C C C 15 Cos �45 700 1e-4 1.5e5 71 b 40�N C W 32�56 15 L �45 700 1e-4 1.5e5 Zero heat flux on southern half 72 b 40�N W C 64�28 15 L �45 700 1e-4 1.5e5 Zero heat flux on western half 73 b 40�N C C C 15 R 25:2�C 700 1e-4 1.5e5 QH=(T*-T) � 35W/m2/K (Haney 1971) 74 b 40�N C C C 15 D exp.73 700 1e-4 1.5e5 Init.: T@equilibrium from exp. 73 75 b 40�N C C C 15 D exp.73 700 1e-4 1.5e5 Init.: T@equilibrium+perturbation 76 b 40�N C C C 15 D exp.73 700 1e-4 1.5e5 Init.: T=4�C 77 f 40�N C C C 15 R 25:2�C 700 1e-4 1.5e5 QH=(T*-T) � 35 W/m2/K (Haney 1971) 78 f 40�N C C C 15 D exp.77 700 1e-4 1.5e5 Init.: T@equilibrium from exp. 77 79 b 40�N C C C 15 R 29:-2�C 700 1e-4 1.5e5 QH=(T*-T) � 10 W/m2/K 80 b 40�N C C C 15 L �45 100 1e-4 1.5e5 MOM Isopycnal mixing 700m2/s 81 b 40�N C C C 15 L �45 0 1e-4 1.5e5 MOM Isop.mixing700m2/s+GM1000m2/s 82 b 40�N C C C 15 L �45 700 1e-4 1.5e5 MOM GM90 1000m2/s 83 b 40�N C C C 15 L �45 700 1e-4 1.5e5 Wind forcing 84 b 40�N C C C 15 L �45 700 1e-4 1.5e5 Seasonal cycle in surf. heat flux Table 2. Variability diagnostics for the experiments detailed in Table 1. exp # Specific. (short) [m2/s] Oscil. Period [year] Oscil. Index [10-3K] KE mean .1kg/m/s2 KE st.dev .1kg/m/s2 MO mean [Sv] MO st.dev [Sv] PHT mean [PW] PHT st.dev [PW] Tsurface mean [�C] Tsurf st.dev [�C] Tbott. mean [�C] Tbott. st.dev [�C] 0 b 40�N 25.4 41.8 .286 .046 8.08 0.90 .254 .025 13.776 .091 3.183 .0012 1 Spherical 23.3 36.8 .254 .038 8.44 0.96 .239 .022 13.443 .084 3.227 .0009 2 Cartesian 25.1 34.4 .286 .045 8.16 0.90 .253 .023 13.740 .087 3.201 .0012 3 f 40�N ~24.9 96.3 .289 .083 7.99 1.85 .259 .051 12.034 .182 3.352 .0049 4 f 20�N ~17.4 94.7 .401 .127 12.21 5.02 .266 .065 11.047 .231 3.516 .0063 5 f 30�N ~20.2 98.5 .345 .111 9.13 2.78 .264 .059 11.592 .209 3.424 .0079 6 f 50�N ~25.9 89.3 .260 .064 7.42 1.43 .257 .037 12.433 .146 3.277 .0086 7 f 60�N 73.0 121.0 .225 .048 7.49 1.40 .257 .026 12.556 .158 3.274 .0037 8 KH= 350 ~25.0 118.4 .337 .079 9.65 2.47 .270 .053 13.833 .211 3.193 .014 9 KH= 500 ~21.2 131.5 .312 .101 8.85 3.00 .271 .064 13.873 .260 3.139 .0093 10 KH=1000 28.8 32.7 .239 .028 9.69 0.67 .250 .022 13.644 .085 3.238 .0015 11 KH=1500 31.1 50.6 .191 .021 11.90 0.78 .243 .021 13.361 .096 3.315 .0023 12 KH=2000 30.7 38.4 .161 .009 12.39 0.58 .234 .011 13.015 .052 3.365 .0016 13 KH=2500 0.0 0.0 .137 .000 12.69 0.02 .223 .000 12.709 .001 3.403 .0000 14 KH=3000 0.0 0.0 .117 .000 12.95 0.00 .213 .000 12.468 .000 3.427 .0000 15 ups.hor.ad ~32.5 62.9 .434 .035 8.90 0.50 .255 .016 13.901 .101 3.194 .0025 16 KV=1e-5 16.9 62.3 .085 .023 3.12 1.17 .264 .064 32.114 .400 3.163 .0099 17 KV=3e-5 16.6 1.8 .150 .001 4.76 0.04 .252 .001 21.291 .014 3.205 .0012 18 KV=3e-4 35.1 84.4 .539 .132 12.96 2.43 .254 .039 9.538 .091 3.171 .0022 19 KV=1e-3 42.7 200.0 1.072 .388 22.07 8.54 .258 .092 6.782 .126 3.264 .0154 20 No Conv ~33.0 107.0 .346 .104 9.50 3.03 .262 .051 12.670 .292 3.225 .0130 21 f40 NoCo ~65.0 125.4 .316 .111 8.48 3.30 .263 .042 10.864 .367 3.342 .0156 22 MOM 26.1 .263 .046 9.26 .282 13.632 3.188 23 PGLslip 37.4 105.0 .294 .062 8.33 1.83 .260 .040 13.213 .211 3.243 .0015 24 PG0 32.3 56.0 1.033 .152 8.70 1.04 .248 .020 13.318 .078 3.094 .0011 25 PG0slip ~45.7 29.6 .587 .047 10.86 1.20 .291 .062 13.167 .074 3.321 .0038 26 PGR0 24.6 29.3 .215 .027 8.55 0.52 .254 .015 14.646 .064 3.167 .0011 27 PGRslip 44.3 85.2 .180 .029 8.65 1.61 .254 .026 14.097 .157 3.264 .0020 28 PGRW 34.0 30.1 .302 .025 7.14 0.65 .255 .015 13.405 .101 3.348 .0003 29 PGS 0.03 34.8 .291 .025 7.03 0.65 .254 .015 3.331 30 PGS 0.3 0.0 0.0 .275 .000 6.81 0.00 .254 .000 3.332 31 AH=1e3 31.6 48.8 .914 .135 8.43 0.97 .248 .020 13.359 .078 3.100 .0011 32 AH=1e4 29.7 43.6 .594 .095 7.85 0.74 .249 .021 13.450 .084 3.128 .0011 33 AH=5e4 28.6 40.4 .393 .071 7.66 0.79 .251 .024 13.652 .094 3.148 .0010 34 AH=1e5 26.9 37.8 .325 .056 7.81 0.84 .252 .025 13.749 .094 3.165 .0011 35 AH=5e5 22.7 35.0 .194 .018 12.61 1.36 .257 .022 13.774 .089 3.277 .0012 36 AH=1e6 20.4 36.1 .150 .015 15.05 2.39 .258 .028 13.943 .121 3.285 .0005 37 AV=1e-4 26.1 .263 .046 9.26 .282 13.632 3.188 38 AV=1e-3 26.0 .263 .045 7.05 .239 13.570 3.188 39 AV=1e-2 26.0 .261 .045 8.33 .260 13.726 3.187 40 AV=5e-2 25.7 .249 .039 7.96 .246 13.904 3.192 41 AV=7e-2 25.8 .245 .038 6.78 .224 13.823 3.195 42 16�14 28.3 25.3 .336 .032 9.94 0.93 .266 .018 13.721 .074 3.267 .0014 43 24�21 26.4 44.1 .311 .058 7.99 1.22 .258 .032 13.857 .112 3.169 .0010 44 48�42 32.8 36.7 .261 .011 9.12 1.10 .250 .017 13.435 .043 3.258 .0011 45 64�56 26.2 39.7 .266 .036 9.26 0.96 .250 .020 13.418 .083 3.215 .0015 46 nz=90(50) 26.3 40.8 .281 .048 8.16 0.97 .253 .026 13.705 .095 3.165 .0009 47 nz=30 26.1 39.7 .282 .047 8.28 0.96 .252 .025 14.073 .093 3.171 .0008 48 nz=20(25) 26.0 38.0 .284 .048 8.04 1.02 .253 .026 14.085 .095 3.179 .0011 49 nz=18 25.9 38.2 .283 .048 8.02 0.99 .253 .026 13.712 .095 3.179 .0011 50 nz=15' 26.1 38.3 .281 .048 8.02 1.01 .254 .026 13.025 .093 3.179 .0010 51 nz=8 11.0 21.8 .315 .004 10.36 0.77 .254 .017 13.894 .014 3.228 .0007 52 nz=8' 14.3 34.8 .294 .012 9.58 0.52 .255 .010 13.007 .013 3.228 .0007 53 nz=8'Upz 33.0 87.0 .436 .117 16.56 4.30 .258 .045 11.399 .153 3.243 .0022 54 nz=4:Upz 43.8 143.0 .745 .174 32.21 6.01 .259 .039 10.324 .156 3.530 .0070 55 nz=4'Upz 40.2 124.0 .654 .162 27.78 5.89 .256 .040 9.494 .111 3.442 .0066 56 nz=2:Upz 83.3 191.0 2.128 .277 25.14 3.12 .270 .034 6.485 .070 3.854 .0041 57 nz=2'Upz 63.7 211.0 2.228 .520 35.89 6.09 .270 .048 5.720 .095 3.785 .0120 58 QH= 5 100.5 7.2 .030 .002 5.90 0.44 .026 001 5.695 .013 3.792 .0006 59 QH= 15 57.6 32.7 .098 .017 7.42 0.96 .083 .011 8.184 .066 3.572 .0013 60 QH= 30 35.4 33.9 .193 .031 7.75 0.82 .168 .018 11.173 .076 3.348 .0013 61 QH= 60 ~18.6 93.0 .363 .082 9.23 2.15 .346 .058 16.122 .194 3.046 .0032 62 QH= 90 ~17.8 182.0 .520 .147 10.88 3.41 .536 .108 20.572 .381 2.782 .0163 63 QH=135 ~16.5 212.0 .747 .183 12.77 3.54 .810 .150 26.386 .470 2.532 .0280 64 LX 50% 19.3 33.8 .193 .041 5.62 1.21 .129 .017 12.736 .120 3.392 .0008 65 LX150% 29.0 51.7 .370 .038 11.90 1.39 .380 .033 14.589 .065 3.026 .0015 66 LX200% 38.1 85.0 .456 .077 14.41 2.01 .501 .065 15.307 .097 2.844 .0030 67 be0.hf 37.5 56.5 .195 .043 11.57 2.99 .251 .050 9.811 .107 3.273 .0013 68 bn0.hf 37.0 87.4 .159 .048 8.52 2.78 .262 .069 9.267 .135 3.484 .0034 69 f40symco ~35.0 248.0 .386 .093 6.83 2.21 .338 .074 13.716 .318 2.886 .0190 70 f40.coshf ~28.0 136.0 .331 .117 8.21 2.51 .326 .067 13.096 .257 3.320 .0142 71 bs0.hf 29.0 56.0 .300 .048 21.03 5.32 .253 .055 11.184 .096 2.885 .0050 72 bw0.hf 28.4 27.5 .167 .005 17.56 1.38 .242 .008 11.940 .005 3.175 .0007 73 25:2�R35 0.0 0.0 .275 .000 15.45 0.00 .222 .000 13.499 .000 3.539 .0000 74 diag. flux 0.0 0.0 .275 .000 15.45 0.00 .222 .000 13.499 .000 3.539 .0000 75 df.pert. 23.0 56.5 .290 .065 12.68 4.40 .224 .042 13.889 .107 3.423 .0014 76 df.Init4�C 23.0 56.6 .289 .064 12.67 4.41 .224 .041 13.890 .107 3.424 .0012 77 f40R25:2� 0.0 0.0 .384 .000 21.73 0.00 .290 .000 13.500 .000 3.688 .0000 78 f40diag.fl. ~19.3 80.5 .363 .055 16.93 5.09 .288 .035 13.725 .129 3.584 .0081 79 29:-2�R10 28.7 33.0 .238 .026 9.52 1.32 .203 .022 13.515 .038 4.671 .0144 80 ISOPYC. ~20.1 .354 .077 12.25 .327 14.293 3.112 81 GM90 15.6 .165 .010 8.26 .242 13.859 3.312 82 GMHOR. 20.8 .136 .020 7.81 .208 13.633 3.212 83 PGL wind 25.1 50.0 .668 .097 7.55 1.74 .253 .040 12.744 .077 3.085 .0011 84 Seas.Cyc. 25.0 49.0 .306 .114 10.52 3.12 .258 .035 12.918 .207 2.879 .0054