ALAIN COLIN DE VERDIERE AND THIERRY HUCK
Laboratoire de Physique des Oceans, Université de Bretagne Occidentale, Brest, France.

submitted to Journal of Physical Oceanography, August 1st 1997.

Corresponding author adress: A. Colin de Verdière
Laboratoire de Physique des Océans, Université de Bretagne Occidentale
6 Av. le Gorgeu, B. P. 809, 29285 Brest cedex, France
Email: acolindv@univ-brest.fr

CONTENT

ABSTRACT

1. Introduction

Interdecadal variability of the temperature of the surface layers of the ocean has been identified very early by Bjerknes (1964) and later by a number of workers, Deser and Blackmon (1993), Kushnir (1994), Hansen and Bezdek (1996). The variability takes the form of large scale, surface intensified anomalies. The EOF analysis of Deser and Blackmon for the North Atlantic reveals the dominance of a dipole clearly intensified near the western boundary roughly situated off Newfoundland. Reverdin et al. (1997) showed that the anomalies are surface intensified and that the salinity signal in the North Atlantic is coherent with the temperature signal. They point out that the region of the slope current is a likely origin for the long period fluctuations. Following Stommel's (1961) idea that different boundary conditions for temperature and salinity could lead to multiple steady states of the thermohaline circulation (THC), Marotzke (1990), Weaver and Sarachik (1991) noted that oceanic GCMs run under mixed boundary conditions (restoring for temperature, flux for salinity) can also exhibit decadal oscillations with clear advective origins. The study of the oscillations was then pursued in simpler contexts with only one active variable forced by surface constant flux by Huang and Chou (1994), Greatbach and Zhang (1995), Cai et al. (1995), Chen and Ghil (1995). Greatbach and Zhang pointed out the strong similarity between the oscillations observed in such square box ocean model forced solely by constant heat flux and those of the fully coupled GFDL ocean-atmosphere model described by Delworth et al. (1993). The oscillations obtained under constant flux were shown to persist in a similar form by Chen and Ghil (1996-CG hereafter) when a simple ocean model was coupled to an Energy Balance Model (EBM) of the atmosphere, removing the assumptions of constant flux at the air-sea interface. What was demonstrated in that study is that at the very low frequencies which are of concern here, the ocean sees almost a constant heat flux with time variations smaller than mean values. A sizable fraction of the constant solar flux at the top of the atmosphere drives the ocean below. Furthermore a number of sensitivity runs allowed them to suggest that the transition from a steady THC to an oscillatory one occurred through a Hopf bifurcation as either the atmospheric turbulent heat diffusivity or the ocean-atmosphere coupling coefficient decreased.
Although the advective origin of the interdecadal oscillations under mixed boundary conditions was early recognized by Weaver and Sarachik (1991), the precise mechanisms under which the oscillations proceed have remained elusive. Winton (1996) demonstrated clearly the three dimensional character of the oscillations by comparison with two dimensional simulations that did not exhibit this type of interdecadal variability. Winton (1996) and Greatbach and Peterson (1996) suggested that key to the oscillatory nature of the THC is the existence of boundary trapped waves that propagate in the Kelvin wave sense around the basin. The Winton's chain of arguments describing the oscillation by the trigger of thermal anomalies through thermal wind currents impinging normal to the coast and propagating in an uninterrupted manner around the basin is difficult to follow and there is no identification of energy sources necessary to sustain the waves against dissipation. Greatbach and Peterson observed that only the western boundary was crucial to the existence of the oscillations and suggested that southward propagating boundary waves perturb the western boundary current that in turn generate perturbations that are advected to the North East corner and play a role in reinitiating the wave propagation. It is not clear to us how the boundary waves excited at the North East corner can make their way back to the western boundary given the absence of stratification along the northern boundary for the boundary waves have velocity that depends on such stratification. Furthermore the explanation seems inconsistent with their demonstration of the absence of sensitivity to the strong relaxation along the eastern boundary. Huck et al. (1997-HCW hereafter) launched a series of experiments under constant flux to show that forcing amplitudes and rotation were conducive to the oscillatory state (stronger amplitudes) while the mixing processes such as convection, horizontal mixing and dissipation were of a damping nature (smaller amplitudes). The presence of the western boundary current was also found to be crucial as shown earlier by Greatbach and Peterson. Runs without the b effect showed that the variation of the Coriolis parameter with latitudes was similarly unimportant as noted earlier by Winton. Although the character of the oscillations is essentially three dimensional, HCW showed that the oscillation can be described using only two active variables: the strength of the overturning and the meridional South-North temperature difference anomaly . Through heat conservation, the rate of change relates to the temperature advection while it is observed in the numerical experiments that itself does not relate to as in Stommel's (1961) model but instead that relates to the existence of a phase lag between overturning and temperature gradient was also observed by Greatbach and Peterson. It is the existence of this quadrature between the "interior forcing" and the western boundary current (which is the dominant contribution to the overturning ) that sustains the oscillation. The origin of this time delay is really three dimensional through either advection or interior westward wave propagation and is reminiscent of the spin up of an ocean basin driven by a wind stress curl (Anderson and Gill 1975), in which the western boundary current builds up and narrows on the time scale it takes for Rossby waves to cross the basin. In the present context b and so Rossby waves have been shown to be unessential to the oscillations. One may expect, however, that given the presence of an underlying mean circulation and stratification, more general potential vorticity (PV) waves should exist in pure f-plane cases to propagate the thermal anomalies between the interior and the western boundary and allow the system to oscillate. While these PV waves are ubiquitous in some of the simulations that we have analyzed, we want to draw the attention in this paper to "the wavemaker" that must be present to sustain the waves in the face of dissipation because the coarse resolution models that we are using are indeed very dissipative! We propose that baroclinic instability of the western boundary current region is such a wavemaker. To demonstrate this proposition a number of planetary geostrophic (PG) simulations are carried out in a situation where the ocean is first put to steady equilibrium under restoring boundary conditions at the surface and then allowed to depart from this initial state under the diagnosed flux kept constant from there on. This procedure generates the familiar western intensification of buoyancy flux that is a key feature of large scale ocean-atmosphere interaction. Once in this configuration of constant flux, classical stability analyses of the numerical solutions are possible and the dominant destabilizing factors identified (Section 2). If our conjecture that baroclinic instability is the dominant process driving the oscillation is correct, then the simplest model exhibiting baroclinic instability at planetary scales is a 3 layer model. On these long interdecadal time scales, the barotropic mode is always in equilibrium with its forcing and it is then unlikely that it may play an important role. It is, in any case, absent from our flat-bottom, buoyancy driven experiments with no bottom friction. Both analytical calculations of the growth rates and numerical simulations under constant flux are carried out in the 3 layer model showing indeed the importance of the first and second baroclinic modes in organizing the instability (Section 3). Finally we continue CG's lines of thought in Section 4 to construct a coupled ocean-atmosphere box model that shows that their important suggestion of a transition from a steady state to oscillatory states via a Hopf bifurcation can be interpreted in a somewhat different way: we propose that the limit cycles of the oscillations arise when the growth rate of the baroclinically unstable western boundary current region exceeds the diffusive time scale resulting from the cumulative damping action of turbulent oceanic-atmospheric diffusivity and infra-red back radiation. On these interdecadal time scales, the atmosphere is restricted to a damping role with its variables enslaved to the active oceanic variables.

2. Evidences for baroclinic instability

a. Critical damping terms
HCW's experiments have identified the parameters that govern interdecadal oscillations, with essentially the amplitude of the surface flux driving the oscillations and all mixing processes but vertical mixing acting as a brake. Figure 1 summarizes some of these results that demonstrate the geostrophic/inviscid character of the oscillations. The reader is referred to HCW's paper for details of the experiments that led to the construction of these figures. Essential parameters of these high resolution (HR) runs are recalled in Table 1. We first define an oscillation index as the basin mean standard deviation of temperature over one period of the oscillation and consider the influence of an imposed surface heat flux zonally uniform and linearly varying in latitude. As the amplitude of the forcing increases, so does the oscillation index (Fig. 1a). If instead the forcing is kept constant but the vertical mixing inceases, the oscillation index again increases (Fig. 1b). Since higher vertical mixing implies a stronger overturning of the mean circulation, both results stress that the more energetic the mean circulation, the larger the amplitude of the oscillations. As shown in HCW, the horizontal diffusion appears as the main damping for controlling the oscillation amplitude (Fig. 1c). For values beyond 2500 m2s-1, no oscillation is observed. From this critical value of the damping, we can infer the actual growth rate of the instability that sustains the oscillations. At the horizontal resolution (160 km) of these sensitivity experiments (HR runs), a time scale Æx2/KH of 120 days emerges. In comparison the vertical mixing plays a very weak role to damp the anomalies: The diffusive time scale Æz2/Kv reaches 3 years for a 100 m depth interval at the pivot value of 10-4 m2 s-1. This very weak direct damping effect of vertical mixing is in fact more than counterbalanced by the increase in the mean THC that follows an increase in Kv. Of course the convective adjustment acting as a very large vertical diffusion where static instability occurs plays a significant damping role on the oscillations as shown in HCW. If the surface forcing is now a relaxation towards an apparent atmospheric temperature, we can again find a critical value for the relaxation coefficient that separates the oscillatory solutions from the steady ones (Fig. 1d). The critical value so obtained corresponds to a time scale of 105 days which agrees well with the one inferred above from variations of the horizontal diffusion. This is a valuable result in view of the analysis of Seager et al. (1995) who showed through modelling of the response of the lower atmosphere to SST anomalies that the sensitivity of heat flux with respect to SST was of the order of 15 W m-2 K-1, significantly less than values currently used in ocean models that come from an assumption of infinite heat capacity for the atmosphere. Such a value translates to a restoring coefficient of order 1 cycle/year for 100 m vertical resolution and the present results suggests that Seager et al.'s inferred values are low enough to allow interdecadal oscillations to occur. To summarize the damping nature of horizontal diffusivity and the active nature of the mean state advection, the Peclet number UL/KH apears a key parameter that controls the strength of the oscillations (the sensitivity of heat flux to SST which is also very important could be measured against advection by a number such as ).

b. Description of the oscillation
We describe here the spontaneous oscillations that arise after the model has been spun up for thousands of years to steady equilibrium under restoring boundary conditions, the surface fluxes diagnosed and kept constant from there on. Having many runs of different resolutions at our disposal, we have decided to illustrate the oscillation with a low resolution run that can be easily reproduced with modest computing equipment (LR run, see Table 1). Of course most of the analysis that follows has been reproduced at both low and high resolution with no essential differences. After an initial transient phase of 10-15 years (described later) the perturbations grow and evolve into an oscillatory state (Fig. 2) that is independent of the initial perturbation that triggers the instability. The thermal anomalies are surface intensified and prominent in the north-western corner where the surface fluxes are the largest. The current anomalies circulate along the contours of the temperature anomalies leaving little net eddy heat flux except near the western boundaries. The motions are geostrophic with little interior divergence. In much the same way as for the mean field, the anomalies of the divergence field are concentrated along the western and northern boundaries, upwelling and downwelling occurring there to connect the surface current anomalies to the deep current anomalies (of the opposite sense since no net barotropic transport is permitted). In figure 2, the oscillation is shown initially in a state of weak thermal structure. On the other hand, the WBC anomaly (and the overturning) is strongly positive, being fed through upwelling along the western boundary. Within about 6 years a major western intensified positive thermal anomaly has developed in response and covers half of the basin. The WBC anomaly is now oriented in the northward direction the mass being supplied to a larger extent by the southern branch of the warm anomaly. The net overturning is weaker and the situation has therefore evolved from a prevailing vertical recirculation phase to a horizontal recirculation phase. In the following nine years, a westward propagation of the warm anomaly is apparent. As the anomaly hits the western boundary, the WBC anomaly disappears and reverses when the southward moving branch of the warm anomaly has reached the western wall. At this point, the vertical recirculation phase has resumed but in the opposite direction with little horizontal interior recirculation. Then the southward WBC creates an offshore interior cold anomaly and the second half of the period of oscillation proceeds similarly.
Characteristic phase diagrams (Fig. 3a and 3b) in the x-t plane shows a partition of the domain between a western third where the temperature oscillations are large and stationary and the interior in which weaker thermal anomalies propagate westward against the mean eastward circulation in the northern part of the basin. In the meridional plane the WBC anomalies present also the character of stationary oscillations intensified in the northern part of the domain (Fig. 3c). This description raises many questions about the choice by the dynamics of the anomalies temporal and spatial scales but we defer a discussion of these to the end of this section and concentrate on the factors that drive and brake the oscillations. It is readily apparent that the drive of the oscillations is an instability in the northern part of the mean western boundary current that builds the vertical recirculation along the western boundary. Consider the initial stages of figure 2, once warm waters have invaded the northern interior through a combination of horizontal advection-diffusion: the anomalous thermohaline circulation is driven uphill against the meridional pressure gradient since anomalous high pressures are associated with the warm anomalies. It is during the phase of propagation of the interior warm anomaly that the THC is defeated and braked. But at the end of this phase, southward current anomalies are present along the western boundary that are seen to amplify. To support this idea that the driving of the oscillations is an instability of the WBC region, we have calculated the terms that govern the growth or decay of the variance of the temperature anomalies (figure 4). The largest positive term appears to be the down gradient eddy meridional temperature fluxes that dominates in the northern third of the western boundary. Since the mean is negative, positive eddy fluxes at the western boundary are at the heart of the existence of the oscillations. Fortunately in these experiments, the regions with positive values (enhancing the anomalies) remain in the same location such that the time-averaged pattern is similar to instantaneous situations (on the contrary under zonally uniform flux, temperature anomalies do travel around the northern half of the basin and the region where the terms are positive do vary a lot along a period: the time averaged pattern is then almost an order of magnitude smaller than the instantaneous snapshots and there is no well defined driving area). The vertical structure of the eddy fluxes (Fig. 5) shows a surface and western intensification. Given that is nearly geostrophic, the occurrence of eddy fluxes relies on variations in the vertical of the perturbations. The comparison of the two zonal sections of figure 3 for temperature and velocity first shows that the zero velocity contours (pressure extrema) are displaced offshore of the temperature extrema. It is this situation of pressure extrema to the east of temperature extrema that allows the meridional velocity to correlate significantly with the temperature anomaly. Because the pressure integrates vertically the temperature field under the hydrostatic constraint, this correlation is therefore possible only if the temperature distribution presents phase lags in the vertical direction. We show in figure 6 the temperature field as a function of depth and time in the unstable regions diagnosed by the positive meridional heat fluxes. A strong vertical phase lag appears between the upper levels dominated by convection and the lower levels that lag by a quarter period. The fact that no such vertical shift is observed in the stable regions of the basin reminds us of the 3D organisation of unstable waves found under quasi-geostrophy. This familiar index of vertical phase lag is indeed necessary to allow downgradient eddy heat fluxes and release of potential energy when the flow is in approximate geostrophic and hydrostatic balance. These observations strongly suggest that the basic driving mechanism of the oscillation is a local baroclinic instability of the western boundary region. An energetic confirmation (not shown) comes also from the nearly perfect phase opposition between the total potential energy and the kinetic energy (whose most part comes from the overturning). This association of vertical phase shift of the temperature distributions and baroclinic instability helps to understand why the oscillations are observed to be damped by the convection scheme adjustment (HCW) since the latter acts to remove all vertical structure and phase lags necessary for the instability.
Given the horizontal boundary layer structure of the mean velocity profile near the western wall, and the complicated distributions of vertical shears necessary to equilibrate the mean buoyancy loss at the surface, the basic state that we must deal with is far more complicated than the idealized zonal flows situation commonly used in baroclinic instability studies. However the signature of the unstable motions of concern here are in a state of near geostrophy and under such dynamics there are not many ways of organizing an unstability. The theoretical stability analysis of a mean state such as ours being a daunting perspective we judge the efficiency of this baroclinic instability process by neglecting entirely the horizontal boundary layer structure of the mean state. We shall gauge the growth rates of the instability against dissipation in the simplest situation that we can envision: a 3 layer model with mean vertical shears only (Section 3).

c.The initial instability
Since we allude to an instability mechanism as a way to sustain the oscillation against dissipation, it is worth looking at the transient phases immediately after the constant buoyancy flux have been switched on. When the previous run (with restoring surface boundary conditions) has been integrated for 3-4000 years, it appears that such a state is stable under a switch to flux conditions when initial external perturbations are absent. A large range of initial perturbations (uniform temperature anomaly, different Fourier modes) have been added to show that this transient stage is eroded over a 10-15 year time scale. The transient state itself is patterned after the initial perturbations and for instance with a uniform surface temperature anomaly, a front appears along the boundary of the convection region, but progressively the anomaly builds up in the north-western quadrant to take the organized structure of figure 2. The transients are then replaced by the previously described oscillatory state which usually reaches steady amplitudes after several oscillation periods. We have checked that the initial growth rates of the perturbations during the first few years (as measured from the rms velocity or rms temperature), is independent on their initial amplitudes and that there are no threshold value that needs to be overcome (the case for linearized intability calculations).
On the other hand we have noticed that if the run with restoring boundary condition is further from equilibrium due to a shorter time integration, there is no need for external perturbations to trigger the oscillations following a switch to flux conditions. After only a few months, a perturbation builds up in the western boundary region where the buoyancy loss (negative heat flux) is largest. This is in the unstable region that was pointed out previously. We may expect that any imbalances in that region will amplify rapidly because of the baroclinic instability mechanism which is strongest where mean vertical shears are largest. To equilibrate heat losses at the surface by horizontal heat advection of nearly geostrophic currents requires precisely the existence of such large mean vertical shears.

d.The structure of the temperature anomalies
The surface intensified temperature anomalies that emerge under flux boundary conditions after a transient phase of several oscillation periods have a well defined spatial and temporal structure that needs to be rationalized. The first obvious comment is that away from the western boundary current region, the interior perturbations are geostrophic and neutrally stable with fluid circulating along the isotherms in agreement with the weakness of the eddy fluxes (Fig. 4). To understand why the anomalies are not passively advected by the dominant eastward mean flow in the northern part of the basin, let's consider a warm anomaly embedded in a mean temperature gradient: west (east) of its center, northward (southward) eddy velocities bring warm (cold) water that propagates the anomaly to the west. This is exactly the classical argument for Rossby waves propagation with the mean meridional temperature gradient taking the role of the effect, the difference being that a mean eastward flow is also associated with the mean temperature gradient to oppose the westward propagation. In some of our experiments (LR runs), the westward propagation dominates while in others, the two effects equilibrate. The interior anomalies that we observe are therefore potential vorticity waves, modes of a surface intensified eastward jet. It can be shown analytically in the inviscid f-plane case of the three layer model developed in the next section that among the two1 eigenvalues, one vanishes and the other is equal to the mean velocity of the middle layer, therefore much less than the mean eastward velocity of the surface layer. In most of the simulations with diagnosed flux, these modes are neutral but in the simulations with zonally uniform imposed flux, these modes can also be unstable in the interior eastward current region. It is not our objective here to calculate these modes in details, only to show how their propagation characteristics help to explain the structure of the observed anomalies. To do so the various terms contributing to the heat equation for the anomalies have been analysed. In the interior, horizontal terms dominate over vertical ones for both advection and diffusion. A zonal section of these terms at a given instant (Fig. 7) reveals that the leading terms governing the evolution of thermal anomalies are mean horizontal advection of anomalies, eddy advection of the mean temperature field and horizontal diffusion. The two advective terms dominate over lateral diffusion and oppose each other, with usually the eddy advection slightly bigger, a situation which in that case conforms with westward propagation. In such instance this means that the interior anomalies could be modeled to zeroth order with: (1) where (positive for westward propagation) would be the eigenvalue associated with the PV modes previously mentioned. Adimensionalizing by the width of the basin and time by , (1) transforms into: (2) with , the Peclet number. If we assume that the waves cross the basin in 10 years, the Peclet number of the LR simulation turns out to be of order of 25, markedly greater than one. Introducing a Fourier solution as in (2) leads to : (3) When the Peclet number is large and and of order one, the two roots of (3) are approximately: The root is the westward propagating wave while the root is an eastward propagating wave that decays to the east as . The observed interior anomalies (western intensification, westward interior propagation) are patterned from these two elementary waves that are forced by lateral boundary conditions. The westward mode is apparent in the interior (figure 3a) while a hint of the eastward propagation of the mode can also be seen near the western boundary. The prediction of the eastward decay scale, a trapping scale for the PV waves turns out to be dimensionally , which varies from one to two grid points in the simulations, to larger values because of the uncertainty in the values of . The equation (1) relates the interior westward propagation of the anomalies to their western intensification over this scale in a manner that mirrors the usual vorticity discussion of Rossby waves type perturbations for oceans at rest. Finally the period of the oscillations is probably the most difficult to rationalize. When the westward propagation is observed, half the period of the oscillation corresponds very crudely to the travel time of the PV waves across the basin. If the size of the basin and the mean state are the two important parameters that in turn determine the eigenvalues of the PV modes then the period would vary roughly as (HCW's work showed that horizontal diffusion does not influence the period of the oscillations). The application of this result remains however somewhat speculative since there are cases where the interior propagation is not so obvious. Furthermore HCW's experiments indicate a scaling for the period closer to . Future work should clarify this point and explore how the structures of the interior PV waves that exist for a given mean state bear upon the period that emerges.

3. Baroclinic instability in simplified analytical and numerical models

Motivated by the previous idea that baroclinic instability of the western boundary current region provides the energy source for the oscillations to be maintained against dissipation, we wish to inquire from a more general theoretical ground the character of that instability and evaluate its strength as gauged by the growth rate. Although baroclinic instability is of central importance to eddy production in mid-latitudes oceans and atmospheres at the scale of the Rossby radius of deformation, the question is really whether it may be active at scales much beyond the Rossby radius. The Phillips's (1954) two layer model of a zonal flow shows that this is not the case because of the existence of a large scale cut-off that suppresses the instability of the largest scales. With such a crude vertical resolution, the barotropic and baroclinic modes are incapable of strong interactions essentially because their time scales become too unequal, the barotropic time scale becoming very much smaller than the baroclinic one for large spatial scales. Models with continuous vertical structure such as Eady's or Charney's do not have such large scale cut-off essentially because higher baroclinic modes are involved in the interaction. Studies of the THC at coarse resolution in PG models neglect entirely relative vorticity an assumption that is justified by the excellent comparison with full PE models (Huck et al. 1996). The consequence is that the barotropic mode becomes entirely diagnostic and hence unable to participate in the baroclinic instability (in fact in our box geometry, buoyancy driven simulations devoid of bottom friction, it is exactly zero). Therefore we are led to think that large scale baroclinic instability is aloowed if interactions between higher baroclinic modes are possible and that models of the THC under constant flux would perhaps behave in a rather different way if the number of layers in the vertical were smaller or larger than two, a point to which we will come back. This question of baroclinic instability at planetary scales has already been adressed theoretically by Colin de Verdière (1986) and Cavallini et al. (1988). In the inviscid case under PG dynamics, Colin de Verdière showed that the presence of b (and horizontal divergence) was crucial to the existence of the instability and that only westward shears were unstable, the growth rate increasing linearly with wavenumber. Friction was added to bound the growth of high wavenumbers. Although this instability is theoretically possible, none of the thermohaline circulation simulations reported so far in the litterature ever mentioned active baroclinic instability. Given that the unstable perturbations are surface intensified, we believe that the reason lies with the presence of the usual strong restoring boundary conditions on surface temperature and/or salinity. When the surface fields are restored on a time scale shorter than the growth rate of the baroclinically unstable perturbations, there is no way that the instability may be triggered. The situation is vastly different however under flux boundary conditions (or mixed boundary conditions) because no such external controls exist to damp the unstable waves. What is then needed at this point to confirm this idea is to find out the favored spatial scales, the growth rates, the form of the perturbations for large scale flows in approximate geostrophic balance and judge the vigor of the instability by comparing with the previous numerical calculations. Although dissipation is small, we are looking for a mechanism that applies for beta-planes as well as f-planes and in the absence of relative vorticity, a process such as friction that allows to break the geostrophic constraint is needed. Given that the instability appears in the simulations in regions such as near the western boundary where friction, although small, is not negligible, leads to study the problem under PG dynamics with Laplacian friction included. Of course the experiments show the instability of a very complicated basic state, that of a western boundary current hugging a wall with vertical as well as horizontal shear-a situation far too complicated to study analytically. Consequently we simplify the picture and consider the stability conditions with respect to vertical shears only.

a. Linear instability calculations As argued previously, the simplest model that can demonstrate the instability is the three layer model. The three layers are of density and thickness with -3 numbering the layers from the top. A rigid lid imposed at the surface and a flat bottom gives: (4) where is the uniform fluid depth. Since the barotropic mode is absent, the condition reduces the problem from 3 to 2 degrees of freedom in the vertical. The pressure being hydrostatic, the geostrophic velocity can be expressed in terms of the layer thickness after use of the two above conditions (k is here a vertical unit vector): (5) Because the two upper layers are supposed to represent the main thermocline, and have been assumed small compared to . The reduced gravity and . Although most studies of baroclinic instability choose basic states which consist of zonal flows, the experiments do not suggest this to be a particularly good choice since the predominantly meridional flows along the western boundaries appear to be potentially the most unstable. Departing from the zonal flow assumptions implies on b-planes to study the instability of forced flows. However when the buoyancy forcing is stationary, the case of interest here, the mass conservation equations for the perturbations do not contain the forcing and hence for each layer the linearized perturbations (with primes) obey: (6) where and represent the velocity and layer thickness of the mean state. A first remark concerns the last term , which is usually not present in baroclinic instability calculations, the mean flow being assumed divergenceless. However suppose that the mean flow is horizontally convergent (as would happen in the neighborhood of solid boundaries for instance), then we can expect exponential growth of the layer at a rate which is just . Take the example of the surface layer in a cooling region to illustrate the effect of that sole term: if the full depth becomes larger than the mean depth then the convergence of mass flux exceeds the constant diapycnal mass loss to the lower layer (the buoyancy forcing) and the lower interface will continue to deepen showing an unstable situation. Although the divergence of the mean flow in the numerical experiments is non zero, an examination of the heat balance has shown that this term is small compared to the first three and therefore unlikely to play a significant role in the instability. It is a local analysis that assumes a constant mean flow that is made in what follows and we neglect entirely mean flow divergence from now on. In the inviscid case (i.e. with geostrophic velocities) it is not difficult to show that a necessary condition for linear perturbations (varying as ) to be unstable is that the quantity changes signs between layers one and two, where is simply the mean PV. If the mean PV gradients are parallel, then instability requires them to be opposite in each layer, while if they rotate from one layer to the next it is always possible to find a direction of a wavenumber vector that allows the quantity to change sign. This necessary condition appears to be a rather weak one because it is fullfilled in many regions which do not exhibit instability in the numerical experiments. To find out in practice where and how strongly baroclinic instability can be, the consideration of friction is essential. Adding a Laplacian dissipation in the horizontal momentum equations makes a correction to the geostrophic perturbation velocities as: (7) The velocity amplitudes above are function of time only and is an horizontal, wavenumber dependent, Ekman number. Knowing for each wavenumber how the perturbation thicknesses relate to velocities allows to cast equations (6) in terms of thickness, from which it is then easy to calculate eigenvalues. To complete the calculation, a mean flow profile must be chosen. The most unstable region appears to be slightly North of the mid basin position of the western boundary region. The surface current has a north eastward direction that veers sharply to the left going deeper, a veering that is entirely consistent with the warm water advection by the WBC to fill the requirement of surface cooling there. The thermocline then deepens to the East and even more to the South due to the large zonal shear that is observed (Fig. 8). With these mean flow velocities typical of the western boundary current region where instabilities grow most rapidly, figure 9a shows that typical growth rates of 1 cycle per year are obtained. The wavenumbers maximizing the growth rate decrease with increasing friction and with the friction coefficient typical of the numerical experiment, one expect growth rates at scales (estimated as quarterwavelength) of the order of 5 to 10 Rossby radius of deformation. For a given friction coefficient, figure 9b shows that the growth rate and the wavenumber maximizing the growth rate increase almost linearly with the upper vertical shear. These simple calculations show that the predicted growth of anomalies in the western boundary current region is rapid compared to the observed period of the oscillations and to the diffusive time scale of the anomaly that was shown to be a major damping factor. They are therefore robust enough to be the unstable modes that are at work in the constant flux experiments.

b. Layered model numerical calculations. To further assess the efficiency of the instability process, a number of exploratory runs have been carried out with the same layered formulations, the diabatic forcing being introduced through mass exchange between the layers. To do so, the mass conservation equations for each layer are advanced in time and the velocities are then computed from the height fields assuming PG dynamics (Rayleigh friction was used for simplicity). We selected rather small forcing amplitudes to prevent surfacing of the layers so that the solutions are quasi-linear and remain close to the regime of the previous analytical calculations. When wind forcing is absent, the solution is internal, so that a two layer case (with a small upper layer) becomes equivalent to a one and a half layer model. We started experimenting with such a model and found that we could never reproduce an oscillatory solution. Although it is possible in principle to cause an instability of a different nature through the divergence of the mean flow (the last term in (6) already discussed), its growth rate O(1cycle/16 years) is too small and does not exceed the critical value imposed by the horizontal mixing. The next step was then to add an intermediate layer, necessary for baroclinic instability to occur (Colin de Verdière 1986). With this two and half layer configuration, we managed to reproduce unsteady behavior with typical variability periods around 25 years. The important point that we want to stress lies with this fundamental difference between one and a half and two and a half layer models regarding the variability under constant flux: while the one and a half layer model never induced any variability after a few hundred years of integration, the two and a half layer model driven by exactly the same forcing within the same geometry exhibited irregular decadal fluctuations of significant amplitudes. These numerical experiments support our proposal that long wave baroclinic instability is the 'driving' for decadal variability in ocean models forced by quasi-steady surface buoyancy fluxes. 4. A coupled ocean-atmosphere box model Given first that the ocean is observed to oscillate on interdecadal time scales under constant heat flux and second that the driving of the oscillation is linked to the instability of the western boundary current region, we propose to go a step further in discussing the implications of this mechanism for coupled ocean-atmosphere models. The coupled model study of Delworth et al. (1993) shows variability on a 50 years time scale whose origin is associated with variations in the intensity of the THC, resulting in western intensified large scale SST anomalies that bear encouraging similarities with oceanic observations. Indeed GP made the case that what Delworth et al. saw in their model was oceanically driven with atmospheric perturbations following the oceanic perturbations generated by a mechanism similar to what is found in coarse ocean only models. Intermediate in complexity, CG's study confirmed that the oscillations persist almost unaltered in the situation of an ocean model coupled to an atmospheric EBM. The objective that is pursued in this section is to use our knowledge of the phase relations between the oceanic meridional temperature gradient and the overturning of the THC, to explore through a coupled ocean-atmosphere box model the conditions under which what has been observed in the ocean model can persist with an atmosphere overhead. Of course the box models have no quantitative ambitions but allow instead to wrap up interactions of large scale feedbacks that may be sought after later in GCM's. Consider then the situation in figure 10, where two atmospheric boxes are coupled to two oceanic boxes. The atmospheric boxes exchange heat externally through incident solar flux and infrared back radiation flux to outer space. We assume that the heat flux between ocean and atmospheric boxes is equal to , i.e. proportional to the difference between ocean and atmosphere temperatures. The characteristic response time of the atmospheric to an oceanic thermal anomaly is , ratio of atmospheric thermal inertia to the coupling coefficient. While a value of 40 W m-2 K-1 for the coupling coefficient is customary in ocean modelling (following Haney 1971), Seager et al. (1995) have suggested that this is too large and that values of 15 W m-2 K-1 or less are more appropriate..With such values and an height of the tropopause , , the response time of the atmosphere turns out to be of the order of one week, so small compared to the characteristic oceanic interdecadal time scales that it is an excellent approximation to assume that the atmosphere is always in equilibrium balance with its fluxes. This assumption is likely to be robust against any known uncertainties in the coupling coefficient , so that for each atmopheric box we write: (8) where , the turbulent heat diffusion coefficient (in W m-2 K-1), parameterizes the turbulent exchange by the large scale atmospheric eddies. The heat conservation equation for the oceanic box is: (9) where is the horizontal area of box , and are respectively the THC overturning strength and the oceanic heat diffusivity, the two agents that strive to homogenize the heat content between the two oceanic boxes. Since our main interest is not in the mean state but in the interdecadal variability, we assume, in the following, a known mean state and concentrate our attention on the meridional temperature differences respectively for the ocean and for the atmosphere. The choice of a time scale (= 1 year) allows to rewrite the equations (8) and (9) for the temperature differences as: (10) and (11) where is the difference of solar flux between the tropical and polar box, the infrared flux has been linearized around a mean state as , is , and is now scaled by . From (10), one obtains immediately: (12) showing that the atmosphere temperature differences depend linearly on oceanic temperature differences, being reduced by a factor that illustrates the dissipative roles of large scale atmospheric turbulence and infrared back radiation. Values of and are indeed significantly smaller than , yielding a reduction of atmospheric temperature differences by about 10 %. The ocean-atmosphere heat flux difference between the boxes can then be expressed entirely in terms of the oceanic variables as: (13) Quite naturally the flux driving the ocean is made of 2 parts, a constant part that is slightly reduced compared to the solar flux at the top of the atmosphere and a dissipative part that summarizes the effect of atmospheric turbulence and infrared back radiation to damp the oceanic anomalies. The time scale involved in this damping is . For thermal anomalies of depth O(500 m), values of and respectively of 2 and 0.6 , a value of 200 years is obtained, showing the weakness of the temperature dissipation that the ocean model sees. Furthermore we recover the point emphasized by CG namely that the first constant term is far greater than the second one justifying the use of a constant flux at very large scales. When is eliminated, a single equation for , the oceanic temperature anomaly is obtained: (14) Where the last term of this equation is: The parameter summarizes the dissipation of the oceanic temperature anomalies through both oceanic and atmospheric eddies and infrared back radiation. Remarquably enough, each of these three effects are observed to have the same order of magnitude1 in the present state of the climate system. Linearizing (14), the sought after thermodynamic equation for the deviation from the temporal mean state (denoted by an upper bar) is: (15) Consistent with the remark of HCW that in the numerical experiments , the term has been neglected. We agree that it would be desirable to have a box model that would correctly represent both the mean and the variable part but unfortunately we have not been able to do this and therefore leave this assumption as a necessary adjustement of the box model by the 3D model results and take the mean as obeying these asumptions. To close the system we need another equation for the dynamics. Originating from Stommel's (1961) study, the box models that have attempted to describe the mean state usually assume a linear diagnostic relation between the overturning and the meridional density gradient. Based on observed behavior of the 3D models, a time lag exist between the overturning and the interior meridional temperature gradient so that the relation between the overturning and the density gradient becomes prognostic (HCW): (16) The second term that has been added on the RHS of (16) models the linear growth rate due to baroclinic instability in the western boundary current region while the last is a saturation amplitude limiting term2 that drives the system back to stability at large amplitudes. The combination of equations (15) and (16) forms a dynamical system in the plane of the two active variables and (we drop the primes for simplicity) that we argue produces the qualitative physics of interdecadal oscillations in the coupled ocean-atmosphere system. We assume that and are given and study this simple system as a function of the dissipation (typically smaller than 1) and of the overall linear growth rate of the instabilities. The last parameter is simply there to parameterize the stabilizing effect of the perturbations at large amplitudes. The PoincarŽ-Bendixon theorem (see for instance Nayfeh and Balachandran, 1995) indicates that a necessary condition for the existence of periodic solutions is that the divergence takes both signs in the plane: So that periodic solutions can exist only if is greater than . Now when the product is smaller than unity, there is no mean state other than zero in these perturbations equations, a desirable situation since we study deviations from the mean. The linear stability properties near the origin can be studied assuming perturbations varying as and we obtain the eigenvalues equation: (17) in which the notation introduces the oscillation frequency of the system when and are zero. When , a condition satisfied in the weakly dissipative, weakly unstable case that we study, it is readily seen from study of the roots of (17) that must be larger than for the origin to be unstable. In this case the roots have also an imaginary part so that the origin is a spiral source. A finite amplitude limit cycle appears and Hopf bifurcation occurs in this parameter range of small dissipation as becomes larger than . As the state system goes away from the origin, the limiting term becomes important and brings the system back toward the origin. Choosing a mean state meridional temperature contrast of 10°C, , and an oceanic turbulent diffusivity of , leads to an overall coefficient of dissipation . The amplitude of the oscillation is plotted on figure 11c as increases from zero while keeping all other parameters constant. The familiar behavior indicates that a genuine Hopf bifurcation occurs at . Once the mean state is chosen, the constant fixes the period. With the above parameters, the limit cycle of the solution is shown in figure 11a and 11b for . Such values of for the box model represent growth rates averaged over the 3D domain in the numerical experiments that will be normally much less than what baroclinic instability can produce locally. It is this low value of the overall growth rate needed to sustain finite amplitude oscillation in the box model that make such prototype equations particularly satisfying ones. They suggest that the instability mechanism of the western boundary current region of oceanic basins should have no difficulty in generating SST anomalies on a time scale short compared to the cumulative dissipation effects by atmospheric-oceanic turbulence and infrared back radiation. On these interdecadal time scales the atmosphere reacts passively, enslaved to the periodicities of the oceanic temperatures (figure 11a). CG conjectured such a Hopf bifurcation from the square root dependency of the THC amplitude in their numerical model when the atmospheric diffusivity and the coupling coefficient were varied. Both the decrease of and lowers the overall dissipation parameter in our box model and therefore increases the THC amplitude reponse because in each case the damping action of the atmosphere is less. What we have added to the picture of CG is that the properties of the limit cycle depend upon the instability of the mean state as measured by the parameter and the dissipation .

5. Concluding remarks

The present analysis of coarse resolution ocean models suggests that western boundary current regions may be baroclinically unstable at scales beyond the Rossby radius in the real ocean and drive interdecadal oscillations. In ocean models driven by constant buoyancy fluxes with sufficiently low subgrid scale diffusivity, phase lags appears (i) between temperature and pressure in the horizontal direction and (ii) between temperatures in the vertical direction. Such phase lags are precisely what is required to produce downgradient eddy heat fluxes under nearly geostrophic dynamics. These fluxes release the mean potential energy in the ouflow regions of the surface western boundary current and may therefore be associated with the wavemaker of the interdecadal oscillations. Excited by this source, surface intensified potential vorticity waves of interdecadal period appear in the more stable interior. With the help of a three layer model, analytical baroclinic instability calculations show that the growth rate of the unstable perturbations in western boundary current regions is compatible with the stability boundaries that delineate the presence of oscillations with respect to the horizontal diffusivity or the restoring time constant (for runs with restoring boundary conditions instead of constant flux). We have advanced the conjecture from the results of an idealized ocean-atmosphere box model that these free oscillations generated through an instability process in the ocean (the high heat capacity fluid) may easily force oscillations in the atmosphere (the low heat capacity fluid). It appears that the present levels of turbulent eddy activities in both fluids (whose effect is to damp the oscillations) are consistent (sufficiently low) with a persistence of the oscillations. Athough we hope that we have described the processes in a convincing way at the large scale, we believe that the decisive step will be to show that similar physics persist when the turbulence at the Rossby radius of deformation is explicitly resolved (Spall 1996). We know that such turbulence has a role in the ocean which goes far beyond what a simple parametrization through a diffusion law can produce. Similarly nonlinear interactions of unstable waves and zonal flows in the atmosphere cannot be excluded as a source of low frequency variability of their own as numerical experiments (James and James 1989, James et al. 1994) or thermally driven experiments in rotating annulus geometries (Fruh and Read 1997) tend to indicate. There are indeed some studies that take such low frequency atmospheric variability as granted to propose that the interdecadal oceanic response amounts to a simple integration of atmospheric white noise (Frankignoul et al. 1997).
Of course there are several other mechanisms that have been put forward to rationalize the presence of interdecadal variability. Weaver et al. (1993) have emphasized that such variability can exist under mixed boundary conditions if the E-P flux has sufficient amplitude. In the context of the present study, this may be understood from the viewpoint that if the E-P amplitudes are large enough, buoyancy anomalies could be generated through an instability process similar to that we have discussed herein for temperature. Because the E-P latitudinal distribution tends to drive a reverse cell, the mean state with the combined thermohaline forcings is weaker than with thermal forcing alone, leading to a weaker baroclinic instability process. Equally important will be to judge the efficiency of the present mechanism when mechanical forcing at the air-sea interface is included. Strong additional feedbacks exist, either local as between the oceanic mixed layer and the wind stress and/or remote as between the western boundary current tranport and the wind stress curl. Such processes have been considered initially by Bjerknes (1964) from analysis of observations and extended more recently by Latif and Barnett (1994) from analysis of coupled GCMs. We feel confident that the robust mechanism that we have identified in our idealized ocean models might play an important role in more realistic contexts (see Capotondi and Holland 1997 for instance).

Acknowledgments This is a contribution of the UMR 6523, partnership between Cnrs, Ifremer and Université de Bretagne Occidentale.

REFERENCES

Anderson, D. L. T., and A. E. Gill, 1975: Spin up of a stratified ocean with applications to upwelling. Deep Sea Res., 22, 583-596.
Bjerknes, J., 1964: North Atlantic air-sea interactions. Adv. Geophys., 10, 1-82.
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.
Capotondi, A., and W. R. Holland, 1997: Decadal variability in an idealized ocean model and its sensitivity to surface boundary conditions. J. Phys. Oceanogr., 27, 1071-2-1093.
Cavallini, F., F. Grisciani, and R. Mosetti, 1988: Bounds on the eigenvalues of the planetary scale baroclinic instability problem. Dyn. Atm. Oceans, 12, 71-80.
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.
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.
Frankignoul, C., P. Muller, and E. Zorita, 1997: A simple model of the decadal response of the ocean to stochastic wind forcing. J. Phys. Oceanogr., in press.
Fruh, W. G., and P. L. Read, 1997: Wave interactions and the transition to chaos of baroclinic waves in a thermally driven rotating annulus. Phil. Trans. R. Soc. Lond., 355, 101-153.
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.
Haney, R. L., 1971: Surface thermal boundary conditions for ocean circulation models. J. Phys. Oceanogr., 16, 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.
Huang, R. X., and R. L. Chou, 1994: Parameter sensitivity of the saline circulation. Climate Dyn., 9, 391-409.
Huck, T., A. Colin de Verdière, and A. J. Weaver, 1997 : Decadal variability of the thermohaline circulation in ocean models. J. Phys. Oceanogr. (submitted).
Huck, T., A. J. Weaver, and A. Colin de Verdière, 1996: The effect of different parametrization and boundary conditions applied to the momentum equation in coarse resolution thermohaline circulation models. J. Phys. Oceanogr., (submitted).
James, I. N., and P. M. James, 1989: Ultra-low-frequency variability in a simple atmospheric circulation model. Nature, 342, 53-55.
James, P. M., K. Fraedrich, and I. N. James, 1994: Wave-zonal-flow interaction and ultra-low frequency variability in a simplified global circulation model. Q. J. R. Met. Soc., 120, 1045-1067.
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, 1994: Causes of decadal climate variability over the North Pacific and North America. Science, 266, 634-637.
Marotzke, J., 1990: Instabilities and multiple equilibria of the thermohaline circulation. Ph.D. thesis dissertation, Institut fur Meereskunde, Christian-Albrechts Universitat, Kiel, 126 pp.
Nayfeh, A. H., and Balachandran, 1995: Applied nonlinear dynamics. Wiley Series in Nonlinear Science. Phillips, N. A., 1954: Energy transformations and meridional circulations associated with simple baroclinic waves in a two level quasigeostrophic model. Tellus, 6, 273-286.
Reverdin, G., D. Cayan, and Y. Kushnir, 1997: Decadal variability of hydrography in the upper northern North Atlantic in 1948-1990. J. Geophys. Res. , 102, 8505-8531.
Seager, R., Y. Kushnir, and M. A.Cane, 1995: On heat flux boundary conditions for ocean models. J. Phys. Oceanogr., 25, 3219-3230.
Spall, M. A., 1996: Dynamics of the Gulf Stream/ Deep Western Boundary Current Crossover. Part II: Low-frequency internal oscillations. J. Phys. Oceanogr., 26, 2169-2182.
Stommel, H., 1961: Thermohaline convection with two stable regimes of flow. Tellus, 13, 224-230.
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., J. Marotzke, P. F. Cummins, and E. S. Sarachik, 1993: Stability and variability of the thermohaline circulation . J. Phys. Oceanogr., 23, 39-60. 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, 289-304.

FIGURE CAPTIONS

Figure 1: Sensitivity of the interdecadal oscillations to: (a) Variations of the amplitude of the meridional distribution of surface heat flux. (b) Variations of the vertical diffusivity coefficient (log-log plot). (c) Variations of the horizontal diffusivity coefficient . Note that no oscillations are found for diffusivity coefficients larger than 2500 m2s-1. (d) Variations of the restoring constant (when the surface temperature is restored to a linear meridional temperature distribution). The presence (´) or absence (o) of the oscillations is indicated along with the oscillation index in the former case. The mean overturning strength is plotted on the vertical axis, since the restoring atmospheric temperatures are changed along with the restoring constant. Note that oscillations disappear for values between 20 and 25 W m-2 K-1. All of these results have been obtained with the HR configuration (Table1) for the solid line, but with a purely geostrophic model (no momentum dissipation but no-slip boundary conditions imposed) described in Huck et al. (1996) for the dashed line. The oscillation index is the basin average of temperature standard deviations over a period.

Figure 2: The anomalies of surface temperature and surface currents during a full cycle of the oscillations (LR run). Basin size is 6000 by 4500 km. Frames are 3 years apart. The overturning maximum (minimum) occur at frames 1,9 (5). The largest temperature and velocity anomalies corespond roughly to half a degree and half a cm/s respectively.

Figure 3: Characteristic diagram in the x-t plane,at a central latitude (a) and (b), in the y-t plane averaged over the western boundary current region (c).

Figure 4: The driving terms of the equation for temperature variance (HR run). Largest values occur for merional fluxes in the North West corner. The mean surface circulation that is superimposed shows that the most unstable region is located where the western boundary current turns eastward. Figure 5: Longitude-depth section of the merional eddy fluxes at a central latitude (LR run). Values must be multiplied by 1e-6 to obtain cm s-1 K units. Figure 6: Characteristic z-t diagram of the temperature field at a point situated in the unstable North-West corner region (LR run). The phase shift of the temperature distribution appears below the mixed layer (the depth axis is horizontal and the contour interval 0.04 C).

Figure 7: The terms of the temperature equations at a given time (year 45 after switch to flux condition, LR run) and at mid-basin latitude as a function of zonal distance across the basin. (a) Vertical terms (b) Horizontal terms

Figure 8: The vertical structure of the mean flow in the unstable regions of the western boundary current for the LR run (u solid, v dashed).

Figure 9: The growth rate (cycle/year) in the 3 layer model as a function of wavenumber (scaled by the inverse Rossby radius of deformation). Parameter values are h1=100m, h2=200 m, H=4200 m, g'3=0.9e-2, g'1=1.8e-2. (a) for different values of the Laplacian friction coefficient (2.6e5, 5.1e5, 7.6e5, 1e6 m2 s-1). Growth rate variations with respect to the friction coefficient is monotonic with small growth rate associated with large friction coefficient at a given wavenumber. (b) for different values of the upper layer meridional velocity (1, 2, 3, 4 cm s-1). All other velocity components are zero. The friction coefficient is set at 1e5 m2 s-1. The growth rate increases monotonically with the shear at a given wavenumber.

Figure 10: The geometry of the four box ocean-atmosphere model.

Figure 11: Results of the ocean-atmosphere four box model: (a) Oceanic (solid) and atmospheric (dashed) temperature anomalies as a function of time (b) The limit cycle in the phase plane of meridional temperature difference and overturning streamfunction for values above critical. (c) Bifurcation diagram of the amplitudes of the oscillations against the growth rate of the baroclinic instability.

TABLE CAPTIONS

TABLE 1
Configurations of the low resolution (LR) and high resolution (HR) runs discussed in the text. The models are described in details in Colin de Verdière (1988) and Huck et al. (1996) respectively. They are based on the planetary geostrophic equations with Laplacian friction closure for a flat-bottomed cartesian beta-plane centered at 40°N and extending from 20°N to 60°N. KH (KV) is the horizontal (vertical) diffusivity, while AH is the horizontal viscosity (compared to which vertical viscosity is safely negligible).
Model Lat.ext. Long.ext. Depth gridpoints gridspacing KH KV AH
(km) (km) (m) lon*lat*ver (km) (m2/s) (m2/s) (m2/s)
LR 4500 6000 4500 12*16*15 375 2000 1e-4 1.0e6
HR 4480 5120 4500 32*28*15 160 700 1e-4 1.5e5

Footnotes:
1 Under the condition of no net horizontal tranport there are only two eigenvalues.
1 Oceanic mesoscale diffusivity observed to be translates to a heat conductivity of 1 for an area A of 2 000 km by 2 000 km.
2 The precise form of this last term is at this point arbitrary and other alternatives exist such as , in which case the dynamical system would be a Van der Pol oscillator in the limit.