Coupling between Wind-Driven Currents and Midlatitude Storm Tracks

A model for the interaction between the midlatitude ocean gyres and the wind stress is formulated for a shallow-water, spherical hemisphere with finite thermocline displacement and the latitudinal dependence of the long Rossby wave speed. The oceanic currents create a temperature front at the midlatitude intergyre boundary that is strongest near the western part of the basin. The intergyre temperature front affects the atmospheric temperature gradient in the storm track region, increasing the eddy transport of heat and the surface westerlies. The delayed adjustment of the gyres to the wind stress causes the westerly maximum to migrate periodically in time with a decadal period. The behavior of the model in a spherical geometry is qualitatively similar to that in a quasigeostrophic setting except that here the coupled oscillation involves oceanic temperature anomalies that circulate around the subpolar gyre, whereas the quasigeostrophic calculations favor the subtropical gyre. Another difference is that here there is a linear relationship between the period of the coupled oscillation and the delay time for the adjustment of ocean gyres to changes in the wind stress. This result departs from the quasigeostrophic result, in which the advection timescale also influences the period of the decadal oscillation.


Introduction
The possibility of midlatitude ocean-atmosphere interactions is exciting, because it brings into play the ocean in a role beyond that of a passive integrator of noisy atmospheric forcing. As such, it allows for the possibility of enhanced predictability in midlatitude weather patterns on decadal timescales.
Although there is some observational evidence of a correlation between decadal fluctuations in atmospheric sea level pressure (SLP) and in sea surface temperatures (SST) anomalies (e.g., Nakamura et al. 1997;Trenberth and Hurrell 1994), the dominant mechanism for the coupling has not been determined. One of the obstacles in analyzing this process, besides the obvious inadequacy of the observational database on decadal timescales, is that different atmospheric general circulation models (AGCM) respond very differently to similarly prescribed SST anomalies [cf. Peng and Whitaker (1999) and the references therein]. The discrepancies in the response are largely due to differences in the models' climatologies and eddy statistics. This is a crucial problem if, as suggested by Peng and Whitaker (1999), the main effect of anomalous SST is to alter the eddy fluxes of heat and vorticity via small changes in the mean vertical shear. Depending on the particular configuration of the climatological storm track with respect to the SST-induced heating anomaly, the anomalous eddy fluxes can reinforce or reduce the perturbation in the mean flow. Because the eddy fluxes are mostly due to transient eddies, this scenario emphasizes the necessity of appropriately resolving or parameterizing the baroclinic processes that maintain the storm track in midlatitudes.
In a previous study, Cessi (2000) has formulated a model that captures a feedback loop between storm tracks and the oceanic currents in which the SST, through their coupling to the atmospheric heat budget, influences the baroclinic eddy activity in the atmosphere and consequently the surface wind stress that drives the wind-driven flows advecting the SSTs. The model couples two simple modules for the ocean and atmosphere; namely, Stommel's model for the ocean gyres and Green's (1970) parameterization for baroclinic eddies for the midlatitude transport of heat and momentum in the atmosphere. The only prescribed forcing in the model is the net absorbed shortwave heat flux at the top of the atmosphere, and although the basic stratifications in the atmosphere and ocean are prescribed, the momentum and heat budgets, based on conservation laws, produce a remarkably realistic climatology-surface westerlies at midlatitudes, with a well-defined storm track forced by surface heat fluxes on the flanks of the intergyre thermal front. Furthermore, the delayed adjustment of the gyres by slowly propagating baroclinic Rossby waves produced a self-sustained oscillation of the atmospheric storm track. VOLUME 14 J O U R N A L O F C L I M A T E Recently, Miller et al. (1998) have shown that the oceanic component of the feedback loop is in place. An increase in the intensity of the westerly winds across the midlatitudes in the late 1970s to early 1980s provided the opportunity for a case study of the oceanic response to a persistent change in the wind stress. Specifically, there is evidence of a spinup of the gyres with a western-intensified thermocline response, together with an SST signal in the Kuroshio-Oyashio extension.
Encouraged by the recent observations, we revisit Cessi's model and inquire into the effects of the idealized model configuration. One unsatisfactory aspect that can be abandoned without sacrificing the simplicity of the model, is the idealized ␤-plane geometry. The planetary scale of the wind-driven ocean circulation suggests that a more realistic spherical geometry be used. On a sphere, the westward phase speed of Rossby waves becomes a function of latitude, an effect that might be important given the quasigeostrophic result that the timescale for Rossby waves to cross the basin sets the period of the oscillation in conjunction with the gyre advection time. In the present study, we have reformulated Cessi's (2000) model by recasting the ocean module in terms of the planetary geostrophic equations. The planetary geostrophic equations, unlike the quasigeostrophic equations can retain the full variation of the Coriolis parameter and are not restricted to small horizontal variations in the thickness of the wind-driven layer.
The plan for the paper is as follows. In section 2, we describe the atmospheric module, comprising vertically integrated, zonally averaged heat (section 2a) and momentum (section 2b) budgets. In section 3, we present the oceanic module based on the planetary geostrophic equations for the momentum budget (sections 3a, 3b) and the thermodynamic budget, based on an advectiondiffusion equation for the SST (section 3c). In section 4, we present the climatology of the model, and in section 5, we describe the variability produced by the coupled model that we contrast with the results obtained from a quasigeostrophic ␤-plane formulation by Cessi (2000). In section 6, a box model is analyzed to clarify the dependence of the oscillation's period on some of the parameters. Finally in section 7, we present a discussion and summarize the results.

The model atmosphere
The strategy of the model is similar to that used in Cessi (2000) except that spherical polar coordinates are used. We make the simplifying assumption that on the timescales of interest, that is, much longer than a month, the atmosphere is in equilibrium with the ocean. Thus, there are two atmospheric diagnostic variables, the surface potential temperature, s , and the surface wind stress, s , that are in equilibrium with the two oceanic prognostic variables, the upper-ocean temperature, T, and the thermocline thickness, h.
In the atmosphere, we consider the zonally averaged, vertically integrated heat and momentum balances. Thus the redistribution of heat and momentum by baroclinic eddies must be parameterized in terms of the zonally averaged quantities. We adopt the parameterizations of Green (1970) and Stone (1972) as detailed in the following.

a. Vertically integrated zonally averaged heat budget
The atmosphere is assumed to adjust instantaneously to the ocean, so that the zonally and vertically integrated heat budget is given by where the angle brackets indicate a zonal average, and is the latitude. Here, Q i is the net absorbed shortwave heat flux at the top of the atmosphere. The Q o is the outgoing longwave radiation that is parameterized by linearizing the ''gray Stefan-Boltzmann law,'' Q o ϭ G() 4 , about a mean value, ⌰ (in kelvins), so that The constants A and B are prescribed, and if the atmosphere were in radiative equilibrium they would determine the mean surface temperature, as well as the difference in temperature between the pole and the equator.
The term F ao is the flux of heat from the atmosphere into the ocean. The zonally averaged air-sea heat flux must be weighed by the fraction of a latitude circle, s, occupied by the ocean. Following Haney (1971), the flux of heat through the ocean's surface is given by the approximation where T s is the SST, and is the bulk heat transfer coefficient. The C pa is the specific heat of the atmosphere at constant pressure, and is the density of the atmosphere, assumed to be a function of height only. The zonally averaged heat flux ͗͘, is parameterized to be down the mean gradient, that is, The form (4) has been adopted by Green (1970) as a plausible representation of the heat flux by midlatitude baroclinic eddies, and has also been derived by Cessi (1998) for the mean meridional circulation in the Tropics. Pavan and Held (1996) summarize the dependence of the eddy diffusivity on the mean heat gradients. Here, for simplicity, we take it to be a function of height only where s and d are constants. The potential temperature profile in the atmosphere is taken to be where the stratification S is a prescribed constant, and the dynamical part of the potential temperature is independent of height. Thus the vertical resolution of the atmosphere is equivalent to that of a two-level model. Last, the vertical structure of the density is taken to be With the above specifications, the heat budget for the atmosphere leads to the following equation for the zonally averaged surface temperature s , The boundary conditions are chosen to ensure that there is conservation of heat in the hemisphere, that is,

. Vertically integrated zonally averaged zonal momentum budget
The longitudinally and vertically integrated zonal momentum balance for the atmosphere in statistical steady state is given by In the above, s is the zonally averaged zonal component of the surface stress, and the left-hand side is the vertically integrated convergence of momentum flux. Therefore, at least when mountain torque can be neglected compared to the viscous surface drag, the surface stress can be estimated through knowledge of the lateral momentum flux. Green (1970), showed that an estimate of the momentum flux can be obtained assuming that it is dominated by baroclinic eddies, ͗u͘ ഠ ͗uЈЈ͘, where primes indicate departures from the zonal average, and making the quasigeostrophic approximation. The momentum transport by the mean overturning circulation, ͗u͗͘͘, is neglected, in accord with the quasigeostrophic approximation. Fortunately, this is not a bad approximation, even in the Tropics, for the vertically averaged momentum flux. To make further progress, the eddy fluxes must be related to zonal mean quantities. Simply parameterizing the momentum flux as downgradient diffusion of the zonally averaged wind is not appropriate because the eddies are observed to accelerate the westerlies in the midlatitudes. However, a local relationship between the potential vorticity and heat eddy fluxes and their respective mean gradients is found for baroclinic eddies growing on large-scale jets (Pavan and Held 1996). Thus, the diffusive closure schemes that apply for conserved quantities when there is a scale separation between the mean and eddy fields, are appropriate for heat and potential vorticity. Green (1970) shows that in the quasigeostrophic approximation, it is possible to relate the flux of momentum to the fluxes of potential vorticity and potential temperature (i.e., heat). For our purpose, it is necessary to apply the ideas of Green (1970), to a spherical geometry, and here we are guided by White (1977) who has generalized Green's model to spherical polar coordinates. An appropriate definition of the quasigeostrophic potential vorticity on the sphere is given by White (1977): where f o ϭ 2⍀ sin o is a typical midlatitude value of the Coriolis parameter, and a is the radius of the earth. In this definition, the planetary vorticity is exactly represented, but the stretching term is approximated by replacing f with f o . This approximation, supplemented by the use of f o in the thermal wind relation for the velocity, o z ⌰ guarantees that the total energy and vorticity are appropriately conserved (Mak 1991). With this definition, the following relationship between eddy fluxes of momentum, potential vorticity, and potential temperature is obtained: Following Green (1970), the eddy fluxes of potential vorticity and of potential temperature are parameterized down the zonal mean gradient, that is, Here, is the same eddy diffusivity used in the atmospheric heat budget, and its expression is given in (5). Substituting (12), (13), and (14) into (9) and using (10), we obtain an equation for the surface stress that contains only zonal mean quantities: Applying the thermal wind relation (11) to the potential temperature profile (6), the zonally averaged zonal wind can be related to the surface temperature and to the surface wind U s through Substituting the vertical profiles, for ͗U͘, ͗͘, , and into (15) and performing the vertical integrations, we obtain a final expression relating the surface stress to the surface temperature profile, The latter constraint can be enforced because in (18) the scale height of the eddy diffusivity, d, is considered to be unknown. In summary, (8) and (18) are the diagnostic equations for the model atmosphere that determine ͗͘ and s once the net absorbed shortwave heat flux, ͗Q i ͘, is specified and once the oceanic surface temperature T s is known. The former is specified in the present model, while the latter is determined by examining the oceanic heat and momentum dynamics.

Model ocean
The ocean model is based on the planetary geostrophic equations in spherical coordinates. This model, unlike the quasigeostrophic equations, can retain the full variation of the Coriolis parameter and is not restricted to small horizontal variations in layer thickness, as is appropriate for the wind-driven ocean circulation on the planetary scale. The planetary geostrophic equations also neglect inertia, and this is a useful approximation for our purposes, since we wish to isolate coupled ocean-atmosphere interactions. Ocean models that retain the nonlinear effects of relative vorticity advection can exhibit intrinsic variability even with prescribed, steady wind stress.
The permanent thermocline is modeled by a single interface separating two fluid layers of uniform but different density. Only the top layer, of thickness h, is set into motion by the wind stress. Below the interface the fluid is assumed to be at rest and isolated from the winddriven circulation. The density jump across z ϭ Ϫh is constant, and gives rise to a reduced gravity gЈ. Although the momentum and mass balance for the ocean model will ultimately be expressed with only one prognostic equation for the thickness of the wind-driven layer, h, the dynamics are more transparent by considering the momentum and mass balances for the Ekman and interior layers separately.

a. The Ekman layer
The layer of thickness h is partitioned into two layers: a top Ekman layer of uniform thickness h e , overlying an interior layer of thickness h Ϫ h e . The dynamics of the Ekman layer is governed by Ϫh ١ · u ϭ w .
e h e e The equations are expressed in spherical polar coordinates, with the longitude, and the latitude. r is the Rayleigh drag coefficient, and the parameter ␣ Ͼ 1, which multiplies the drag coefficient in the meridional momentum equation increases the dissipation in the zonal direction. This anisotropy in the frictional parameterization allows us to adequately resolve the western boundary layer without excessively damping interior perturbations. The choice of a stronger drag coefficient in the zonal direction is also consistent with the numerical simulations of Haidvogel and Keffer (1984), which showed that the variations of the Coriolis parameter with latitude makes eddy diffusivities more effective in the zonal direction. The other parameters are gЈ, the reduced gravity, and a, the radius of the earth. Substituting the horizontal velocity from (21) and (22) into the continuity equation (23) we obtain the expression for the Ekman pumping w e , in terms of the oceanic prognostic variable h and the atmospheric diagnostic variable s .

b. Below the Ekman layer
The dynamics of the thermocline layer, of depth h Ϫ h e , is governed by We have assumed that the bottom of the thermocline, z ϭ Ϫh, is a material surface and thus that there is no exchange of density with the water below the thermocline.
A single prognostic equation in h is obtained combining the continuity equations for the Ekman and interior layers to eliminate the Ekman pumping velocity w e . This equation is most simply expressed in terms of the vertically averaged velocities, (U, V), defined as The vertically averaged velocity is entirely determined by h and s and is given by It only requires knowledge of the wind stress s and it is independent of the Ekman layer thickness h e . We impose the mass-conserving boundary condition that U · n vanishes at the solid walls.

c. Thermodynamics
Although uniform density layers are a useful idealization for the dynamics of the wind-driven ocean, thermodynamics are cumbersome to represent. We thus allow the temperature T to vary horizontally within the layer, without accounting for the horizontal pressure gradients and associated velocity vertical shears that should accompany lateral temperature gradients that are not compensated by salinity.
Thus our oceanic heat balance treats T as a passive scalar, vertically uniform throughout the layer, which is advected by the wind-driven velocity field, (30). Following the treatment in Young (1994) the vertically integrated upper-ocean heat content hT is governed by On the left-hand side, the second term is the divergence of the heat flux carried by the vertically averaged geostrophic plus Ekman currents U, and the third term parameterizes the transport by mesoscale eddies as a diffusive flux down the temperature gradient, ١T. The right-hand side of (32) is the flux of heat through the top of the ocean, given by (3), and C w is the heat capacity of water. We assume that there is no exchange of heat through the bottom of the layer, at z ϭ Ϫh. The oceanic heat budget (32) is guaranteed to conserve heat when supplemented by the insulating boundary conditions A h ١ h T · n ϭ 0 on the solid walls.
With this formulation, the oceanic temperature T is not completely passive, because it affects the atmospheric potential temperature through the air-sea heat fluxes. Thus, through the atmospheric (8) and oceanic (32) heat budgets and the wind stress balance (18), T indirectly influences the velocity field (U, V) by which it is advected. The oceanic velocity is given by (31), and with the specification of the incoming shortwave radiation the model is closed.

Climatology
In this section we describe the climatology of a typical solution, obtained for the set of parameters listed in Table 1.
The only forcing for the coupled system is the incoming net shortwave radiation, ͗Q i ͘ that is constructed by fitting a simple polynomial to the data of Stephens et al. (1981), that is, and is plotted in Fig. 1 (solid curve). Also shown in Fig. 1 is the outgoing longwave radiation, ͗Q o ͘, (dashed curve) that requires knowledge of the atmospheric temperature, and here the time-averaged solution is used. At low latitudes the incoming radiation exceeds the outgoing radiation, while at high latitudes, the opposite is true. To maintain this equilibrium, the wind-driven gyres in the ocean and the parameterized baroclinic eddies in the atmosphere carry poleward the heat absorbed at low latitudes. This increases the surface temperature at high latitudes where the excess heat is radiated back to space. The time-and zonally averaged surface temperature profile is plotted in Fig. 2 (solid curve), along with the temperature profile that would be obtained for a radiative equilibrium solution, (dashed-dotted curve), The temperature at the equator is reduced from the radiative equilibrium temperature by about 10ЊC and warmed at the pole by about 30ЊC, as a consequence of the northward transport of heat by the ocean and the atmosphere.
The time-averaged surface wind stress is plotted in   At low latitudes there are easterlies with a maximum surface wind speed of 4 m s Ϫ1 corresponding to a surface stress of 1.1 ϫ 10 Ϫ1 N m Ϫ2 . For reference, the surface stress that would be obtained if there were no ocean is plotted with the dashed curve. The difference between the two profiles is due to the thermal feedback of the wind-driven ocean heat transport on the driving winds. Clearly the largest effect of the ocean on the atmosphere is in the region of the westerlies' maximum. As detailed in the following, the heat transport by the ocean is northward at all latitudes, and tends on average to reduce the pole-to-equator temperature gradient. However, the reduction in the gradient in the surface temperature field does not happen at all latitudes. The confluent currents in the ocean produce a strong thermal front at the boundary between the subtropical and subpolar gyres that locally intensifies the atmospheric temperature gradient and strengthens the surface westerlies over the corresponding latitude band (cf. the dotted curve in Fig. 2, which is the zonally averaged ocean temperature). Because the surface westerlies are driven by eddy fluxes of momentum, the acceleration of surface flow implies a local intensification of the atmospheric eddy activity by the oceanic flow.
In Fig. 4, the transport stream function for the ocean is contoured. The ocean circulation consists of three counter-rotating western-intensified gyres (there is a weak cyclonic gyre in the Tropics). The maximum transport is 54 Sv (1 Sv ϭ 10 6 m 3 s Ϫ1 ) in the subpolar gyre, and 41 Sv in the subtropical gyre. The meridional transport, driven by the curl of the wind stress in the interior, is returned in swift western boundary currents that are confluent between the subtropical and subpolar gyres. This results in a stretching deformation field that concentrates the isotherms along the intergyre boundary. In Fig. 5, the ocean temperature is contoured. There is a strong temperature front near 55ЊN in the region between the subpolar and the subtropical gyres. In Fig. 6, the heat flux through the ocean surface is contoured, with shading indicating regions of heat flux into the atmosphere. The heat flux field has a dipole structure that straddles the temperature front, with a maximum heat loss to the atmosphere of approximately 450 W m Ϫ2 south of the front and a maximum heat flux into the ocean of 200 W m Ϫ2 north of the front. This distribution of surface heat flux tends to erode the oceanic temperature front. The front is nevertheless maintained by the vigorous western boundary currents that advect cold water southward in the subpolar gyre and warm water northward in the subtropical gyre. In the atmosphere, the converse is true-the surface heat flux pat- 3. The time-averaged surface stress for run 5, (solid curve) and for reference the surface stress that would be obtained if there were no ocean (dashed curve). Although the northward heat transport by the ocean tends on average to reduce the pole-to-equator temperature gradient and the zonally averaged surface winds, the ocean currents produce a strong thermal front at the intergyre boundary that increases the atmospheric temperature gradient and strengthens the atmospheric jet over the corresponding latitude band. tern tends to strengthen the baroclinicity, which is then balanced by the baroclinic eddies that erode the temperature gradient. Thus, the heat fluxes through the airsea interface are responsible for maintaining the increased eddy activity in the atmospheric storm track.
The impact of the air-sea heat flux is best understood by considering the heat budget for the whole system. The time-and zonally averaged heat budget for the ocean is given by

‫ץ‬
ao where the left-hand side is the divergence of the zonally averaged ocean north-ward heat transport. The global heat budget can be obtained by eliminating the surface flux in the atmospheric heat budget (8), multiplying by the area element for a latitude band, 2a 2 cos d, and integrating in latitude: Note that although the factor s appears explicitly only in the ONHT term, its presence is implicit in the ANHT term since s appropriately weighs the air-sea heat flux that determines the zonally averaged surface temperature and therefore atmospheric flow that transports the heat (see also section 6). The zonally averaged nature of the atmospheric module cannot capture the zonal variations of the storm tracks, but their influence on the zonally averaged heat transport are accounted for. In any event, as Thompson and Wallace (2000) show, the dominant mode of atmospheric variability has an essentially zonal structure.
The first term on the right-hand side is the atmospheric northward heat transport, ANHT, and the second term is the ocean northward heat transport, ONHT. The left-hand side is the integrated residual from the equator to between the incoming shortwave radiation and the outgoing longwave radiation. Figure 7 shows the northward heat transports as a function of the sine of the latitude. The oceanic heat transport, (dotted curve) is northward everywhere, with two prominent local maxima at the latitudes of largest mass transport in the subtropical and subpolar gyres. In the subtropical gyre, ONHT peaks at 1.5 ϫ 10 15 W and in the subpolar gyre it peaks at 0.5 ϫ 10 15 W. The atmospheric northward heat transport (dashed curve) peaks at 3.2 ϫ 10 15 W near the latitude of the strong thermal front over the ocean. This maximum in ANHT is reduced and shifted northward from its position for the no-ocean case (dashed-dotted curve). This shift toward the latitude where the ONHT is minimum is due to the reduction in ANHT over the latitude band in the subtropical gyre where ONHT is largest: the atmosphere must transport more where the ocean transports less.
Although many features of the observed climatology are reproduced in our idealized model, others are poorly represented. The ocean temperature in the subpolar gyre is in general too cold. This is probably due to the absence of convection, of sea ice, and of a thermohaline component to the ocean circulation in our model. As a result, the temperature gradient at the intergyre boundary is overestimated, and the resulting air-sea heat flux-   es are closer to winter time values than to the annual mean.

Variability
The climate variability captured by the model consists of a periodic oscillation involving both the atmosphere and the ocean. In the atmosphere, the divergence of the meridional transport of heat and momentum by baroclinic eddies oscillates in response to the slowly evolving sea surface temperature field. The surface wind stress, balanced by the divergence of the momentum transport, also oscillates, and provides a time-dependent forcing to the ocean circulation.
In the ocean, long baroclinic Rossby waves are continually generated near the eastern wall, and propagate westward into the interior, in a continuous attempt to bring the flow into Sverdrup balance with the changing wind stress. The meridional currents induced by the waves perturb the intergyre temperature front and produce SST anomalies that then alter the atmospheric temperature. Therefore the intergyre boundary is the latitudinal band where fluctuations in the oceanic currents can most affect the atmosphere. Moreover, because the SST gradients are strongest in the western part of the basin where the western boundary currents are confluent, the feedback to the atmosphere is delayed by the time it takes for the waves to cross the basin. This delayed feedback prevents dissipative processes in the ocean and in the atmosphere from completely damping out the oscillation.
The typical behavior is illustrated with run 5 from Table 2. The time-averaged fields and heat balances for this run were presented in section 4. In Fig. 8, the anomalous surface wind stress in the latitude band of the intergyre boundary is contoured as a function of time through two cycles each lasting 18.4 yr. Positive and negative wind stress anomalies form near 55ЊN and drift poleward, under the influence of the evolving SST field. The result is an atmospheric storm track that periodi-cally migrates northward before weakening and reintensifying to the south.
The zonally averaged ocean temperature as a function of time is contoured in Fig. 9. The shaded contours are proportional to the north-south temperature gradient. During one period, the temperature front first intensifies near 52ЊN, drifts northward approximately 3Њ of latitude, and then weakens before reintensifying to the south. The perturbation of the intergyre temperature front produces temperature anomalies that are then advected around both gyres. The evolution of the temperature anomalies can be seen in Fig. 10 anomalies greater than 0.2ЊC (shaded) superimposed on the streamfunction field. Alternating cold (dashed contours) and warm (solid contours) temperature anomalies are formed at the intergyre boundary and are then advected by the currents. The strongest anomalies are found in the subpolar gyre, where they propagate eastward across the basin, and are then advected poleward toward the basin boundary, where they are stretched out by the northern boundary current.
In the subtropical gyre, the anomalies are much weaker and dissipate quickly, contrary to the result obtained in the quasigeostrophic calculations (Cessi 2000) where the stronger anomalies circulate in the subtropical gyre. In Fig. 11a we plot the area-weighed variance of the zonally averaged temperature anomalies for the subtropical and subpolar gyres as a function of the control parameter gЈ, and in Fig. 11b we plot the ratio of the two. The ratio of the area-integrated anomalies is close to unity, although the area occupied by the subtropical gyre is over 4 times larger than that of the subpolar gyre. This implies that for all values of gЈ, the local temperature anomalies are stronger in the subpolar gyre. For the particular case of gЈ ϭ 0.022 m 2 s Ϫ1 , the stronger subpolar anomalies can readily be seen from Figs. 9 and 10. Thus, because of the difference in the basin width between high and low latitudes on a sphere, anomalies of comparable area and strength contribute to a different degree to the zonally averaged temperature in the subtropical and subpolar gyres, with the anomalies in the subtropical gyre having the smaller relative contribution. The relaxation to a zonally averaged atmosphere tends to homogenize oceanic temperature zonally resulting in more dilution of an anomaly of given area in the subtropical gyre than in subpolar gyre. This geometrical discrimination is most effective for anomalies generated very close to the intergyre boundary, and cannot be captured by the quasigeostrophic approximation.
Remarkably, Fig. 11b shows that the ratio of the areaintegrated temperature variance in the subtropical versus subpolar gyres drops as gЈ is decreased. Typically, increasing the delay in a delayed differential equations tends to destabilize solutions that are close to neutral (see section 6). Based on this result, one would expect that oscillators in the subtropical gyre that have shorter delays than those in the subpolar gyre might be destabilized by increasing the delay time (decreasing gЈ). Instead, the result that the variance of the zonally averaged temperature rise and fall in unison in both gyres as the control parameter gЈ is varied, indicates that a single unstable delayed oscillator acting at the latitude of the intergyre front is responsible for producing the variance in the temperature field.

a. Heat balance
In the time-dependent zonally averaged heat balance an additional heat storage rate term for the ocean must be added to (36), Figure 12 shows a time series of each term in (37) evaluated at 56ЊN, the latitude where the anomalous ocean northward heat transport is greatest. Two cycles of the oscillation are shown. The maximum anomalous ONHT is more than 3 times larger than the maximum ANHT, but a large part of this transport is stored in the ocean and returned during the next phase of the oscillation. The residual between the transport and storage of heat by the ocean is mostly balanced by the atmospheric heat transport. The radiative imbalance is small throughout the oscillation. The peak in the oceanic heat transport lags behind the peak atmospheric transport by slightly more than a quarter period. In fact, at all lati-tudes in the subpolar gyre, the phase lag for the peak in ONHT is between /2 and . This phase lag is consistent with that obtained with a linear-delayed ordinary differential equation with two feedback terms, one with and one without delay. In such a model, unstable oscillatory solutions must have the delayed feedback term (ONHT) lag the nondelayed feedback term (ANHT) by more than a quarter period, but less than a half period (see section 6).
The ocean heat storage rate can be decomposed into a part due to isopycnal heaving and a part due to isopycnal temperature changes, that is, We have approximated each term on the right-hand side by taking the time average, denoted by an overline, of the undifferentiated field. In this way the roles of the layer thickness changes versus the temperature changes on an isopycnal surface can be isolated. Figure 13 shows a time series of the heat storage term decomposed into the above two terms. As would be expected from a geostrophic scaling the heaving term is significantly smaller than the term due to an isopycnal change in temperature. Hence the fluctuations in the thermocline displacement do not contribute significantly to the time- 12. Time series of anomalies in the ocean northward heat transport (ONHT dash-dotted), the atmospheric northward heat transport (ANHT dashed), the ocean heat storage rate (solid), and the radiative imbalance (dotted), at 56ЊN for run 5. dependent thermodynamic balance. Since the model formulation does not allow a heat flux through the bottom of the dynamically active layer, the large spatial variations in the layer thickness provide a geographically varying oceanic heat capacity.

b. Timescales
The oceanic response to a change in the wind stress, and thus the thermal feedback on the atmosphere, is delayed by the adjustment time for the wind-driven circulation. This delay is equal to the time for a baroclinic Rossby wave to propagate across the basin. In Fig. 14, the position of a wave front originating at the eastern wall is plotted at successive times. The basin crossing time as a function of latitude can be inferred by the intersection of the wave front with the western wall. Variations of the Coriolis parameter and in the mean thermocline depth make the speed of the Rossby waves a function of position. For long Rossby waves, the variation in the thickness h, enters the phase speed only in the Rossby radius of deformation. The contribution from the dependence of the potential vorticity gradient on the layer thickness is exactly canceled by the Doppler shift caused by the fluid velocity that is proportional to the gradient in layer thickness. The phase speed of the long Rossby waves is then simply 2⍀ cos gЈh c ϭ Ϫ .
x 2 a (2⍀ sin) For the basin geometry used in our study, the crossing time, T ϫ , as a function of the latitude is given by where, s ϵ W/360, is the fraction of the latitude circle occupied by the ocean basin. We have conducted two series of numerical experiment in order to explore how the period of the oscillation varies with a change in the timescale T ϫ . In one series of experiments we varied the reduced gravity gЈ, and in the other we varied the width, W, of the ocean basin. Table 2 summarizes the numerical experiments. In Fig. 15, the period of the oscillation is plotted as a function of the crossing time, T ϫ , evaluated at the latitude of the intergyre temperature front. The circles indicate the runs in which the width of the basin was changed and the squares indicate the runs in which gЈ was varied. For the series of runs in which gЈ was varied, there exists a linear relationship between the period of the oscillation and the crossing time. For the series of runs in which the basin width was varied the relationship is also nearly linear except perhaps for a hint of curvature in the runs with the narrowest basins.  Fig. 12, decomposed into a term due to isopycnal heaving (dashed-dotted line) and one due to isopycnal changes of temperature (dashed line).
FIG. 14. Plot of a Rossby wave front at successive times. The front is initially aligned with the eastern wall. The travel time for a Rossby wave to cross the basin at various latitudes can be deduced from the intersection of the front with the western wall.
The two points at the ends of the dashed line in Fig.  15 produced only damped oscillations, they correspond to Runs 4 and 13 in Table 2.
The linear relationship between the period and the delay time differs from that obtained by Cessi (2000) with a quasigeostrophic formulation for the ocean. There, the period of the oscillation was found to scale as the geometric mean of the advective timescale and of the delay. Here, the temperature anomalies advected by the subpolar gyre are diffused along the northern boundary before any fraction can be reinjected at the intergyre boundary. Thus, the advective timescale does not enter in the ocean-atmosphere feedback, which only operates at the intergyre boundary.
In Fig. 15, the slope of the linear relationship between the delay and the period differs depending on whether the delay is varied by changing the width of the basin or by changing the reduced gravity gЈ. We find that increasing the delay by increasing the speed of the Rossby waves through gЈ produces a larger increase in the period than a comparable increase in the delay obtained by making the basin wider. The difference is due to the fact that increasing the width of the basin not only increases the delay, but it also increases the strength of the thermal feedback on the wind stress by increasing the coupling between the ocean and the atmosphere through the parameter s in (8). A simple delayed or-dinary differential equation (see section 6) shows that for a fixed delay, an increase in the strength of the delayed feedback produces a longer period, consistently with eraged deviations from climatology.

A delayed oscillator box model
We illustrate some properties of delayed oscillators by considering a simple delayed ordinary differential equation We motivate the equation with a two-box model for the ocean atmosphere heat budget. Our purpose is not to obtain a quantitative model for the coupled system but simply to illustrate some of the qualitative properties of delayed oscillators. The following heuristic derivation provides a guidance to interpret the different terms in the equation, and to suggest how the parameters vary as the fraction of a latitude circle occupied by the ocean, s, is changed. In particular, we will discuss how the stability, the period of the oscillation, and the phase lag between the oceanic and atmospheric heat transport depend on the parameters. We divide the hemisphere into two regions: one is to the south of the latitude separating the subtropical gyre from the subpolar gyre and the other is to the north of it (Fig. 16). We assume that the changes in the radiative balance can be neglected compared to the changes in the heat transports by the ocean and atmosphere. We then consider the heat budget integrated over each box, and linearized around the steady balance, As in the full coupled model, we assume that the atmospheric heat storage is negligible for the decadal timescales of interest, so that we have the following two slave relations for the atmospheric heat budgets: Taking the difference between (44) and (45) and using the slave relations (47) and (48) to eliminate the airsea heat fluxes F ao1 and F ao2 , we get To close the system, we need to express ANHTЈ and ONHTЈ in terms of the ocean temperature difference ⌬TЈ.
To obtain ANHTЈ, we assume that the atmospheric heat transport due to baroclinic eddies is proportional to the atmospheric temperature difference By combining Eqs. (46), (47), (48) and (52), we can express the atmospheric temperature difference in terms of the oceanic temperature difference The oceanic heat flux anomalies are due to winddriven gyre anomalies acting on the mean ocean temperature gradient. The gyre anomalies in turn are forced by the anomalous wind stress and are therefore proportional to the wind stress at a time t Ϫ t 0 , where t 0 is the the averaged spinup time for the ocean gyre, and is a fraction of the Rossby wave crossing time, T ϫ , defined in (40). Since both the atmospheric transport anomalies and the wind stress anomalies are due to the same physical process, namely baroclinic eddies in the atmosphere acting on the atmospheric temperature gradient, the wind stress is also proportional to the anomalous atmospheric temperature difference. However, the gyre anomalies lag the wind stress anomalies by the gyre spinup time, t 0 . Using (53), we can relate the oceanic heat transport to ⌬TЈ time lagged by the delay t 0 . Finally, the contribution of the oceanic heat transport on the integrated heat balance must be weighed by the fraction of a latitude circle occupied by the ocean basin. This introduces an extra s dependence in ONHTЈ This expression should be compared with (36), where only the s that weighs the contribution of the oceanic heat flux on the zonally integrated heat budget appears explicitly, the other s dependence is implicit in the mass transport term Vh.
With these assumptions, we obtain a single delayed differential equation for ⌬TЈ, the temperature difference between the ocean boxes, Last, we introduce the the following quantities to simplify the notation With the simplified notation, the box model is expressed as a single delayed ordinary differential equation the time for long Rossby waves to cross the basin. It is inversely proportional to the basin width.
To explore the behavior of Eq. (58), we seek solutions in the form . The characteristic equation for the (iϩ)t e frequency and the growth rate , are A graphical solution of t 0 and t 0 as functions of ␣t 0 and ␤t 0 is given in Fig. 17. The contours in the right panel show the growth rate. The thick solid line is the neutral curve separating stable (dashed lines) from unstable solutions (solid lines). For a fixed delay, t 0 , decreasing ␣ or increasing ␤ tends to destabilize the mode. Also, the neutral curve lies above the line ␣ ϭ ␤, and asymptotes to this line as ␣ → ϱ. Hence, stable solutions that are close to the neutral curve will tend to become destabilized if the delay is increased while keeping ␣ and ␤ fixed. The period in multiples of the delay t 0 is contoured on the right panel of Fig. 17. For a fixed t 0 , increasing the strength of the delayed feedback coefficient, ␤ leads to oscillations with shorter periods. Since ␤ is an increasing function of s, the simple box model suggests that increasing s with a fixed delay time should lead to shorter periods. This is consistent with the full model where increasing the width of the basin (i.e., increasing s) for a fixed delay time leads to a shorter period.
One can also deduce the phase lag of the ocean northward heat transport with respect to the atmospheric northward heat transport from the right panel in Fig.  17. From Eq. (56) the phase lag between ONHTЈ and ANHTЈ is given by t 0 phase lag ϭ 2.
(61) period Hence, the contours in the right panel of Fig. 17 give the phase lag, with the lag given by 2 divided by the contour value. The neutral stability curve (thick solid line in panel b) intersects the /2 phase lag contour when ␣ ϭ 0. Hence, the phase lag for unstable parameter values with ␣ Ͼ 0 is always greater than /2 and tends toward, but never exceeds as ␣ and ␤ tend to infinity.
In summary, we have shown how the stability, the period of the oscillation and the phase lag between the delayed and the nondelayed feedback term depends on some parameters. Delayed oscillators that are close to the neutral stability curve will in general become unstable when the delay time is increased. We have shown that the phase lag between the delayed feedback term (ONHTЈ) and the nondelayed feedback term (ANHTЈ) must be between a quarter period and a half period for growing oscillations. We have also shown that the period of the oscillation decreases when the delayed feedback coefficient is increased, and that this coefficient is a monotonically increasing function of the coupling parameter s.

Summary and discussion
In this study we have extended the model introduced by Cessi (2000) for the interaction between the wind and the wind-driven component of the oceanic heat transport. As in the previous formulation, highly simplified and parameterized modules for the ocean and atmosphere are coupled through the global thermal balance. The present model replaces the quasigostrophic formulation with planetary geostrophy in a spherical hemisphere. In the latter formulation the delay time for the adjustment of the oceanic heat transport to a change in the wind stress is a function of latitude. Despite this difference, the delayed negative feedback mechanism VOLUME 14 J O U R N A L O F C L I M A T E identified by Cessi (2000) is still effective at producing a self-sustained oscillation.
In both models, the intergyre boundary is the site at which the ocean-atmosphere coupling is most effective and is where the time-dependent fluctuations in both oceanic and atmospheric variables are localized. This is the region where the time-mean temperature gradients are maximum in both the atmosphere and the ocean, and where the westerlies peak. In this sense, the modeled feedback is between the separated western boundary currents and the midlatitude storm track, through exchanges of anomalous heat fluxes at the air-sea interface.
Although the mechanism for generating the oscillation is the same in the quasi-and planetary geostrophic formulations, two important qualitative differences have been identified. One difference is that in the spherical geometry the temperature anomalies advected into the subpolar gyre are stronger than those circulating in the subtropical gyre, which is the opposite of the quasigeostrophic calculations. The preferred migration of temperature anomalies into the subpolar gyre is due to a geometrical effect on the sphere and it resembles the observations by Sutton and Allen (1997) of the slow northeastward progression of North Atlantic SST. Our calculations indicate that the ultimate fate of the SST generated at the western end of the intergyre boundary is to be advected around the subpolar gyre, suggesting their possible implication in the modulation of deepwater formation in the North Atlantic (Curry et al. 1998).
Another significant difference is that with the planetary geostrophic formulation there is a linear relationship between the period of the coupled oscillation and the delay time for the ocean gyres to adjust to a change in the wind stress. This departs from the quasigeostrophic result where the advection timescale also influences the period of the oscillation. Because of the preferred advection in the subpolar gyre, SST anomalies are stirred against the northern solid boundary where they are rapidly diffused before they can be reinjected in the formation site at the intergyre boundary. Thus advection around the gyre is not implicated in the coupled feedback mechanism. Remarkably, in the more complicated spherical geometry the period has a simpler dependence on the control parameters.
Both advection and Rossby wave propagation have been proposed as possible mechanisms for producing oscillatory behavior in coupled ocean-atmosphere general circulation models. For example, Grotzner and Latif (1998), proposed a delayed oceanic feedback mechanism for the generation of coupled ocean-atmosphere oscillatory modes identified in an extended-range integration with a coupled ocean-atmosphere GCM. However, given the complexity of fully coupled GCMs, the authors have been unable to identify whether advection or Rossby wave adjustment is responsible for the phase switching mechanism they propose. Here we show that oscillatory behavior can be obtained solely because of delayed negative feedback due to Rossby wave propagation.
In this study we have increased the realism of the formulation of the coupled model introduced by Cessi (2000) while retaining much of the simplicity of the original formulation. As such it is an additional contribution to a hierarchy of models of increasing realism necessary for understanding midlatitude interaction between the wind and the heat transport of wind-driven currents. The fact that other sources of variability are absent in the present model made it possible to clearly identify the mechanism responsible for producing the variability. Future work should study how additional sources of variability affect the characteristics of the oscillation and guide us in attempting to detect it in natural data or GCM simulations.