## Abstract

A simple dynamic model is proposed to illustrate the multidecadal oscillation of the Atlantic Ocean thermohaline circulation. The proposed oscillation relies on alternating actions of positive and negative feedbacks, which are operated by a slow adjustment of the ocean circulation and the associated time delay in the advective flux response to a change in meridional density gradient. The key element of the oscillation is the time delay, which is conceptually related to the basin-crossing time of long Rossby waves in the high-latitude North Atlantic. For a sufficiently long time delay, the solution becomes unstable in some regions of model parameter space and oscillates with a period of approximately 2 times the delay time.

## 1. Introduction

The objective of this note is to propose a simple dynamic model for the multidecadal oscillation of the Atlantic Ocean thermohaline circulation (THC), which is commonly manifested in general circulation model simulations (e.g., Delworth et al. 1993; Dai et al. 2005; Knight et al. 2005). By definition, a self-sustained oscillation requires alternating actions of positive and negative feedbacks. Here, it is proposed that the required feedbacks for the multidecadal oscillation of the THC are provided by a slow adjustment of the ocean circulation and the associated time-delayed advective flux in response to a change in meridional density gradient. A four-box model is used to illustrate the proposed oscillation.

## 2. The four-box model

The North Atlantic Ocean is simplified to four boxes with two-layer structures in the high and low latitudes, as shown in Fig. 1. The meridional density gradient is always positive; thus the volume transport must be northward in the upper layer and southward in the lower layer. Conservation of mass dictates that the volume transport at the middepth is downward in the high latitudes and upward in the low latitudes. Therefore, volume integration of the density conservation equation for each box yields

where *ρ*_{1}, *ρ*_{2}, *ρ*_{3}, and *ρ*_{4} are densities of the upper low-latitude box, upper high-latitude box, lower high-latitude box, and lower low-latitude box, respectively; *V* is the volume transport (per unit volume) in response to the meridional density gradient; *q* is density flux into the upper high-latitude box (or out of the upper low-latitude box); *r* is a damping coefficient; *k _{υ}* is a vertical diffusion coefficient;

*H*is the model ocean depth divided by 2; and

*F*

_{1},

*F*

_{2},

*F*

_{3}, and

*F*

_{4}represent other forcing terms such as horizontal diffusion and convective mixing. Note that a cubic form of dissipation is used instead of a linear form to prevent infinite growth of linearly unstable solutions and model drifts as in Suarez and Schopf (1988) for the ENSO oscillator.

Now, a separate equation for *V* is required to solve Eqs. (1)–(4). Using the geostrophic balance and hydrostatic relation, the zonal baroclinic (i.e., upper layer minus lower layer) velocity in the midlatitude corresponding to the north–south density gradient can be written as

where *g* is the gravitational acceleration, *ρ _{o}* is a reference density,

*f*is the planetary vorticity in the midlatitudes,

_{o}*L*is the meridional length of the model domain divided by 2,

_{y}*ρ*

_{14}= (

*ρ*

_{1}+

*ρ*

_{4})/2, and

*ρ*

_{23}= (

*ρ*

_{2}+

*ρ*

_{3})/2. As discussed in Killworth (1985), the baroclinic meridional motions are established after the adjustment time, which depends on the basin-crossing time of long baroclinic Rossby waves. Therefore, in this study, the meridional volume transport (per unit volume) is represented by using Eq. (5) with a time delay:

where *α* is a constant on the order of 1 and *δ* is the time delay, which is approximately the basin-crossing time of long baroclinic Rossby waves.

Scaling time by 2*f _{o}L_{y}*

^{2}/(

*ĝH*),

*ρ*by [

*ĝH*/(2

*rf*

_{o}L_{y}^{2})]

^{1/2}, and

*q*by {

*r*

^{−1}[

*ĝH*/(2

*f*

_{o}L_{y}^{2})]

^{3}}

^{1/2}, the following nondimensional equations can be derived:

where *ĝ* = *g*Δ*ρ*/*ρ _{o}* (Δ

*ρ*is the scale of north–south density difference) and the nondimensional vertical diffusivity

*k*= 2

_{o}*k*

_{υ}f_{o}L_{y}^{2}/(

*ĝH*

^{3}).

## 3. The one-equation model

To gain some insights into the behavior of the nonlinear system of Eqs. (7)–(10), a simplified one-equation model is derived and evaluated here. Specifically, it is assumed that *ρ*_{2}, *ρ*_{3}, and *ρ*_{4} are close enough to each other to be represented as a single variable, *ρ**_{2} = (*ρ*_{2} + *ρ*_{3} + *ρ*_{4})/3. This is not an unreasonable assumption because the convective process maintains *ρ*_{3} as close to *ρ*_{2}, and the advective flux from box 3 to box 4 also keeps *ρ*_{4} close to *ρ*_{3} for a sufficiently small value of *k _{o}*. Then, the density budget equations [Eqs. (7)–(10)] can be simplified to a single equation:

where *ρ̂* = *ρ**_{2} − *ρ*_{1}, and the nonlinear dissipation term is slightly modified from its original form. If *δ* = 0, then this equation has a positive steady-state solution. By neglecting the nonlinear dissipation term for the sake of illustration, the positive steady solution *ρ̂ _{o}* is

The linear stability of *ρ̂ _{o}* in Eq. (11) can be studied by replacing

*ρ̂*(

*t*) in Eq. (11) with the sum of the stationary solution

*ρ̂*and the perturbation

_{o}*ρ̂*′(

*t*). By retaining only the linear terms, the perturbation equation can be written as

Seeking solutions of the form *ρ̂*′ ∝ exp(*σt*) with *σ* = *σ _{r}* +

*iσ*, the following equations can be derived:

_{i}If *σ _{r}* is positive, then the first term on the rhs of Eq. (14) is always larger than the second term; thus Eq. (14) cannot be satisfied. Therefore,

*σ*is negative regardless of

_{r}*σ*. This means that the simplified density budget equation [Eq. (11)] is always stable around

_{i}*ρ̂*. However, it is shown in section 4 that the nonlinear system of Eqs. (7)–(10) is unstable for a sufficiently large value of

_{o}*δ*. This suggests that the density budgets in the lower layers (boxes 3 and 4) must be fully incorporated to resolve the multidecadal THC oscillation.

## 4. Numerical solutions

The behavior of the nonlinear Eqs. (7)–(10) is explored numerically using a fourth-order Runge–Kutta scheme. When the upper layer is heavier than the layer below, the convective mixing is achieved by completely mixing the two layers. Horizontal diffusion is turned off for simplicity because its linear damping effect on the meridional density gradient does not add much value. Figure 2 shows three model solutions for *δ* = 0, 7, and 20, when *α*, *q*, and *k _{o}* are set to 2.0, 0.1, and 0.2, respectively. The dashed lines are the statistical equilibrium values. For small values of

*δ*, the solution achieves a stable state. For a larger value of

*δ*, on the other hand, the solution oscillates with a period of approximately two times the delay time (hereinafter referred to as a delayed advective oscillator).

Obviously, at issue is why the solution oscillates when the advective delay is sufficiently long. To answer this question, the equation for the density gradient between the high-latitude boxes and low-latitude boxes is diagnosed here:

Figure 3 shows (*ρ*_{23} − *ρ*_{14})(*t*) and (*ρ*_{23} − *ρ*_{14})(*t* − *δ*) in the top part of the figure. The storage term on the lhs of Eq. (16) and the advective density gradient flux, surface density flux gradient, and dissipation terms on the rhs of Eq. (16) are shown in the bottom part of the figure. In this case, *α*, *q*, and *k _{o}* are set to 2.0, 0.1, and 0.2, respectively, and the time delay

*δ*is set to 20.

Table 1 summarizes the temporal evolutions of the meridional density gradient, the advective density gradient flux, and the sign of advective feedback during one cycle of the oscillation between the five points labeled (A), (B), (C), (D), and (E) as indicated in Fig. 3. As illustrated in Table 1, the delayed advective oscillator is maintained by alternating actions of amplification (i.e., positive feedback) and stabilization (i.e., negative feedback) through the delayed advective density gradient flux. In one cycle of the oscillation, there are two periods of amplification separated by two periods of stabilization.

The first amplification occurs when the meridional density gradient increases above the equilibrium point because the advective density gradient flux is smaller than the surface density flux gradient [(A)–(B)]. This is followed by the first stabilization period, during which the meridional density gradient swings back from its maximum point toward the equilibrium point because the advective density gradient flux is larger than the surface density flux gradient [(B)–(C)]. Then, the second amplification occurs, during which the meridional density gradient decreases below the equilibrium point because the advective density gradient flux is larger than the surface density flux gradient [(C)–(D)]. During the second stabilization period, the meridional density gradient swings back from its minimum point toward the equilibrium point because the advective density gradient flux is smaller than the surface density flux gradient [(D)–(E)]. This cycle then repeats. In summary, the meridional density gradient and the magnitude of advective density gradient flux anomaly are negatively (positively) correlated during the positive (negative) feedback periods as indicated in Table 1.

To have an idea about the actual time scale of the delayed advective oscillator, let us consider typical parameter values for the North Atlantic: *L _{x}* = 5 × 10

^{6}m,

*L*= 2.0 × 10

_{y}^{6}m,

*ĝ*= 10

^{−2}m s

^{−1},

*f*= 10

_{o}^{−4}s

^{−1},

*β*= 2 × 10

^{−11}m

^{−1}s

^{−1}, and

*H*= 2000 m. According to previous general circulation model simulations, the density fluctuation associated with the multidecadal oscillation of THC is mainly in the high-latitude North Atlantic around 40°–65°N (e.g., Delworth et al. 1993) and propagates slowly to the west (e.g., Dijkstra et al. 2006). Therefore, the long baroclinic Rossby wave speed

*c*(=

*g*′

*Hβ*/(2

*f*

_{o}^{2}), where

*g*′ is the reduced gravity) is computed separately based on typical parameter values for the high-latitude North Atlantic around 50°N:

*g*′ = 5 × 10

^{−3}m s

^{−1}and

*f*= 1.1 × 10

_{o}^{−4}s

^{−1}. Then, the basin-crossing time can be estimate by

*L*/

_{x}*c*, which is approximately 20 yr in this case. Assuming that the delay time is on the order of the basin-crossing time, the period of the delayed advective oscillator is about 40 yr, which is 2 times the delay time. Using the time scale 2

*f*

_{o}L_{y}^{2}/(

*ĝH*) = 1.3 yr, the nondimensional delay time

*δ*is approximately 15.

## 5. Effects of external forcing

Now, the responses of delayed advective oscillator to low-frequency external forcing patterns, such as freshwater flux into the high-latitude North Atlantic, and to high-frequency external forcing patterns, such as weather noise related to the North Atlantic Oscillation (NAO), are explored here. Three low-frequency forcing experiments are performed—that is, 1) density flux out of the high-latitude North Atlantic (i.e., freshening or warming), 2) complete shut down of the THC, and 3) density flux into the high-latitude North Atlantic (i.e., cooling or net evaporation). For each experiment, the following form of density flux is used only in the high latitudes (i.e., box 2):

where *t** = *t*/*δ*, and *q _{o}* is set to −

*q*/2, −2

*q*, and

*q*/2 for experiments 1, 2, and 3, respectively. Note that the surface flux in the low latitude is kept constant in these experiments.

Figure 4a shows that the THC may slow down as a result of the density flux out of the high latitudes (i.e., freshening or warming), but it swings back to the equilibrium solution once the external forcing is removed. An interesting feature is that the amplitude of the delayed advective oscillator is reduced substantially, and its recovery is extremely slow. As shown in Fig. 4b, the THC completely shuts down as the surface flux in the high-latitude box is reduced to match that in the low-latitude box. The THC swings back to the equilibrium solution after the external forcing is removed. However, the delayed advective oscillator is completely disrupted. Interestingly, the external density flux out of the high-latitude North Atlantic (i.e., cooling or net evaporation) has only a minor effect on the THC strength, as shown in Fig. 4c. However, the delayed advective oscillation is weakened (but not as much as in Fig. 4a) and slowly recovers once the external forcing is removed.

Next, the effect of high-frequency forcing on the delayed advective oscillator is explored. It is widely believed that the high-frequency portion of the NAO originates from weather noise. Nevertheless, the NAO has a coherent spatial structure with a dipolelike meridional pattern of the sea level pressure (Hurrell 1995). Because of this coherent spatial pattern, if the high-latitude is cooled (warmed), then the midlatitude is warmed (cooled) during a positive (negative) phase of the NAO. Therefore, the high-frequency forcing of the NAO is represented here as a random noise surface flux with antisymmetric meridional pattern; that is, the sign of random forcing is opposite in the two latitude boxes but with the same amplitude. Note that the random noise forcing does not produce a net surface flux into or out of the system. The amplitude of the random forcing is set to *q*/2.

Figure 5 shows the model solutions with (Fig. 5a) and without (Fig. 5b) the random forcing, and with both the low- and high-frequency forcing (Fig. 5c) for *α* = 1.2, *q* = 0.1, *k _{o}* = 0.2, and

*δ*= 20. In the last case, the low-frequency forcing is added only in the high latitude (i.e., box 2) using Eq. (17) with

*q*= −1.5

_{o}*q*. As shown in Fig. 5b, without the weather noise, the oscillation is damped out for the given parameter values. Interestingly, if the weather noise is introduced, the delayed advective oscillator with the period of ∼2

*δ*can sustain its amplitude of up to about 35% of the mean. It is also interesting to note that, under the weather noise forcing, the amplitude of the oscillation fluctuates at a very low frequency, which amounts to the multicentennial time scale using realistic parameter values for the North Atlantic. Figure 5c shows that, even when the external forcing is large enough to nearly shut down the THC, the weather noise can invigorate the delayed advective oscillator once the external forcing is removed.

In summary, the THC is remarkably stable because it always swings back to its original state once external forcing is removed [this conclusion may not be valid if both the temperature and salinity are considered, as in Stommel (1961)]. The delayed advective oscillator is, on the other hand, very fragile. If external forcing is large enough, then it can virtually wipe out the delayed advective oscillator. The recovery of the delayed advective oscillator is extremely sluggish, suggesting that the growth rate of the delayed advective oscillator is very small. However, these are characteristics of the delayed advective oscillator in its pure form. If the NAO-like weather noise forcing is added to the system, then the behavior of the delayed advective oscillator is drastically changed. In particular, relatively large amplitude of the weather noise (50% of the mean is used in the experiment) can sustain an active delayed advective oscillation of an otherwise stable system. The THC can still shut down if external forcing is large enough. However, the weather noise can quickly invigorate the delayed advective oscillator once the external forcing is removed. Finally, the weather noise can also produces a very low frequency fluctuation of the delayed advective oscillation at the multicentennial time scale.

An important question is why the delayed advective oscillator is excited by the NAO-like weather noise. The simple stochastic climate model of Hasselmann (1976) provides a plausible explanation for this question. It is well known after Hasselmann (1976) that a random noise atmospheric forcing produces a red-noise spectrum of ocean temperature via ocean memory. If this theory is applied to the four-box model, it means that a random surface forcing can produce large amplitude signals in the meridional density gradient field (*ρ*_{23} − *ρ*_{14}) at low frequencies, including at the frequency of delayed advective oscillator *ω* ∼ 0.5*δ*^{−1}. Therefore, the delayed advective oscillation can be excited and maintained by the weather noise even if it is subject to a damped oscillation, as shown in Fig. 5a.

## 6. Linear stability analysis

To better understand how the four parameters *α*, *q*, *k _{o}*, and

*δ*influence the delayed advective oscillator, a linear stability analysis of the nonlinear system of Eqs. (7)–(10) is performed. First, the stationary solutions with

*δ*= 0 are obtained by numerically integrating Eqs. (7)–(10). Replacing the solutions with the sum of stationary solution and perturbation and retaining only the linear terms can derive the perturbation equations (not shown). Seeking solutions of the form

*ρ*′

*=*

_{k}*c*exp(

_{k}*σt*) with

*k*= 1, 2, 3, and 4 and

*σ*=

*σ*+

_{r}*iσ*, a matrix equation, 𝗔 · 𝗖 = 0, can be derived. The determinant of the matrix 𝗔 must vanish for nontrivial eigenfunctions to exist; this yields an equation for the calculation of the complex eigenvalue

_{i}*σ*for chosen values of

*α*,

*δ*,

*k*, and

_{o}*q*. Since the matrix 𝗔 contains the eigenvalue

*σ*and its exponential form, exp(−

*σδ*), an iterative Muller’s method is used to obtain the eigenvalue

*σ*.

Figure 6 shows the neutral curves on the *α*–*δ* plane for *q* = 0.1 and *k _{o}* = 0.2 (Fig. 6a), on the

*k*–

_{o}*δ*plane for

*α*= 2.0 and

*q*= 0.1 (Fig. 6b), and on the

*q*–

*δ*plane for

*α*= 2.0 and

*k*= 0.2 (Fig. 6c). For a given value of

_{o}*δ*, increasing

*α*and decreasing

*k*destabilize the system. These results are not surprising because

_{o}*α*is proportional to the meridional volume transport, which provides the positive feedback required to maintain the delayed advective oscillator and thus serves as a growth rate, and

*k*serves as a damping rate. However, it is important to note that the neutral curve on the

_{o}*q*–

*δ*plane is not monotonic. For a given value of

*δ*, the system is unstable only when

*q*is within a certain range. This strongly supports an idea that the delayed advective oscillator exists via a delicate balance of various terms and also nicely explains why the delayed advective oscillator (in its pure form) is so fragile under the effects of external forcing, as illustrated in Fig. 4.

## 7. Summary and discussion

Perhaps, the four-box model presented here is one of the simplest dynamic models for the multidecadal oscillation of the THC, which is commonly manifested in general circulation model simulations. Despite the overly simplified nature of the model, this minimal complexity model describes effectively the mechanism of the delayed advective oscillator, which appears to be an important stepping-stone toward our understanding of the THC and its multidecadal oscillation.

The key element of the delayed advective oscillator is the time delay in the advective flux response to a change in meridional density gradient. This time delay, which is conceptually related to the basin-crossing time of long baroclinic Rossby waves at the high-latitude North Atlantic, allows alternating actions of positive and negative advective feedbacks and thus gives rise to a self-sustained oscillation. An important signature (or indicator) of the delayed advective oscillator is that the meridional density gradient anomaly and the magnitude of meridional advective density gradient flux anomaly are negatively (positively) correlated during the positive (negative) feedback periods.

Related to the fundamental issue of decadal predictability of the THC, an apparently important and practical question is whether the multidecadal oscillation of the THC is self-sustained (e.g., Greatbatch and Zhang 1995) or damped, in which case the NAO-like weather noise may sustain the oscillation (e.g., Capotondi and Holland 1997; Eden and Greatbatch 2003). The linear stability analysis of the delayed advective oscillator may help us to better understand under what conditions (or regions of the model parameter space for *α*, *k _{o}*,

*q*, and

*δ*) the THC switches between an unstable state and a damped regime.

The delayed advective oscillator can be compared to the delayed action oscillator for ENSO (Suarez and Schopf 1988) because the key element in both cases is the oceanic Rossby wave transit effects. However, there is an important distinction between the two oscillators. In the delayed action oscillator, the equatorial oceanic Rossby wave transit effects provide a negative delayed feedback to an otherwise linearly amplifying system. In the delayed advective oscillator, the advective density gradient flux, which is delayed by the oceanic Rossby wave transit effects, provides both the negative and positive feedbacks for the oscillation.

It is important to point out that, for a given positive (negative) anomaly of meridional density gradient, if there is no time delay then the meridional salt advection has a destabilizing effect whereas the meridional heat transport has a stabilizing effect. Griffies and Tziperman (1995) showed that a phase lag between the salt and heat advection could support a damped oscillation. This mechanism, which is obviously missing in our four-box model, can be readily explored by expanding the density equations [Eqs. (7)–(10)] to two sets of equations for salinity and temperature equations.

## Acknowledgments

We thank R. Greatbatch, Y.-G. Park, and two anonymous reviewers for insightful comments. This work was supported by a grant from the National Oceanic and Atmospheric Administration’s Climate Program Office and by National Science Foundation Grant ATM-0850897. The findings and conclusions in this report are those of the authors and do not necessarily represent the views of the funding agency.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* Dr. Sang-Ki Lee, NOAA/AOML, 4301 Rickenbacker Causeway, Miami, FL 33149. Email: sang-ki.lee@noaa.gov