We provide a dynamical systems framework to understand the Atlantic Multidecadal Oscillation and show that this framework is in many ways similar to that of the El Niño/Southern Oscillation. A so-called minimal primitive equation model is used to represent the Atlantic Ocean circulation. Within this minimal model, we identify a normal mode of multidecadal variability that can destabilize the background climate state through a Hopf bifurcation. Next, we argue that noise is setting the amplitude of the sea surface temperature variability associated with this normal mode. The results provide support that a stochastic Hopf bifurcation is involved in the multidecadal variability as observed in the North Atlantic.
Several pronounced large-scale patterns of variability are known in the present climate system. On the interannual time scale, the El Niño/Southern Oscillation (ENSO) phenomenon provides a dominant pattern in sea surface temperature (SST) variability, which is localized in the equatorial Pacific (Philander 1990). The Atlantic Multidecadal Oscillation (AMO), which is associated with a basin-wide temperature anomaly in the North Atlantic, is the clearest phenomenon on the multidecadal time scale (Enfield et al. 2001). Of course, many ENSO cycles have been measured over the past decades, while the AMO is recorded over one cycle at most. The basic theory of both the ENSO and the AMO should explain (i) the physics of the SST pattern and its propagation (if any), (ii) the physics of the dominant time scale of variability, and (iii) the processes controlling the amplitude of the SST pattern.
For ENSO, such a basic theory exists (Neelin et al. 1998). An important aspect in the development of the theory was the availability of a so-called minimal model of the coupled equatorial ocean–atmosphere system, the Zebiak–Cane (ZC) model (Zebiak & Cane 1987). The solutions of the ZC model have been extensively analysed with the ocean–atmosphere coupling strength μ (the amount of wind stress per SST anomaly) as an important control parameter (Neelin et al. 1998). In the ZC model, the tropical Pacific annual mean climate state can become unstable when μ crosses a critical value μc. When μ>μc, specific time-dependent perturbations grow in time leading to oscillatory behaviour on an interannual time scale.
In terms of dynamical systems theory (Dijkstra 2005), a Hopf bifurcation occurs at μ=μc in the ZC model. The simplest dynamical system (Guckenheimer & Holmes 1990) exhibiting such a bifurcation is the system(1.1a)(1.1b)having two degrees of freedom (x, y). In polar coordinates (r, θ) with x=r cos θ and y=r sin θ, equations (1.1a) and (1.1b) can be written as(1.2a)(1.2b)For μ<μc, there is only one steady state, i.e. r=0 (or x=y=0). For μ>μc, however, there are two solutions of the steady equation (1.2a), i.e. r=0 and . The latter corresponds through (1.2b) to a periodic orbit with angular frequency ω and period 2π/ω. Hence, at the Hopf bifurcation (μ=μc), periodic behaviour with a frequency ω is spontaneously generated through an instability of the trivial solution x=y=0.
In the ENSO theory, when the ZC model is discretized under annual mean forcing, a large-dimensional dynamical system of the form(1.3)appears, where the state vector X consists of the dependent quantities in the model (e.g. SST, oceanic and atmospheric velocities) at each grid point and f contains the tendencies of all these state variables. In the same way as for the simple system (1.1a) and (1.1b), Hopf bifurcations are found by considering the stability of the annual mean Pacific climate state in the ZC model. Putting and , we find by linearization around that (the spatial pattern of the eigenmode) is determined by(1.4)where is the Jacobian matrix of f at . In this case, a Hopf bifurcation occurs when a complex conjugate pair of eigenvalues crosses the imaginary axis as μ crosses μc. If the associated eigenvector is indicated by , then the periodic orbit at the Hopf bifurcation has the form(1.5)and Φ(t) defines a propagating pattern with a period 2π/σi. Note that the period is internally determined by processes in the system and not externally imposed.
As is a solution to (1.4), it is called a normal mode and, in the case of the discretized ZC model, it is usually referred to as the ENSO mode. By following this mode to smaller values of μ, it is found (Jin & Neelin 1993) that it splits up into two modes (in a so-called mode merger or mode splitter). One of these modes (an SST mode) is related to tendencies in the SST equation in the ZC model and the other mode (an equatorial ocean basin mode) is related to equatorial wave dynamics. The pattern of SST in the ENSO mode at μc is inherited from the SST mode while the interannual time scale is inherited from the equatorial basin mode (Jin & Neelin 1993).
While, in the ZC model, sustained ENSO-type oscillations are found when μ>μc, there is no interannual time-scale oscillatory behaviour when μ<μc as the annual mean Pacific climate state is stable. However, when noise (e.g. representing atmospheric weather) is applied in the ZC model, interannual oscillations are found for μ<μc. The simplest system exhibiting qualitatively the same behaviour is the stochastic Itô extension of the simple system (1.1a) and (1.1b), i.e. the dynamical system(1.6a)(1.6b)where λ is the amplitude of the additive noise and Wt is a Wiener process with increment dWt. The expectation value E[Rt], where , resulting from the stochastic integration of the system (1.6a) and (1.6b) is shown in figure 1 for several values of λ; the deterministic case is shown for λ=0. Clearly, there is a response for values μ<μc, which increases with increasing noise level λ.
The effect of noise on the variability in the ZC model has been systematically studied by Roulston & Neelin (2000) and the results are qualitatively similar to those in figure 1. For values μ<μc, white noise in the wind stress over the equatorial Pacific is able to excite the ENSO mode to substantial amplitude, while for values μ>μc there are sustained oscillations that are not much affected by the noise. In both the cases, the spatial pattern and the time scale of propagation associated with the interannual variability do not depend on the precise noise characteristics as both are coupled to the ENSO mode. Stochastic noise in the wind stress mainly affects the amplitude of the pattern (Roulston & Neelin 2000).
Having introduced this stochastic dynamical systems view of ENSO, we now turn to the issue of whether a similar framework applies to the AMO. In fact, the main aim of this paper is to show that such a framework can be developed and that the results are very promising. In §2, we shortly present the minimal model of the AMO and, in §3, we show that a Hopf bifurcation occurs in this model leading to multidecadal periodic oscillations. The origin of the AMO mode giving rise to this variability and the physics of its time scale and pattern is also provided in §3. In §4, the effect of additive noise is presented and we conclude with a summary and discussion of the dynamical systems framework of the AMO in §5.
2. The minimal primitive equation model
A minimal model of the AMO was formulated by Greatbatch & Zhang (1995) and Chen & Ghil (1996) and consists of flow in an idealized three-dimensional Northern Hemispheric sector model forced only by a prescribed heat flux. We therefore consider ocean flows in a model domain on the sphere bounded by the longitudes ϕw=286° (74° W) and ϕe=350° (10° W) and by the latitudes θs=10° N and θn=74° N; the ocean basin has a constant depth H. The flows in this domain are forced by a restoring heat flux Qrest (in Wm−2) given by(2.1)where λT (in Wm−2 K−1) is a constant surface heat exchange coefficient. The heat flux Qrest is proportional to the temperature difference between the ocean temperature T* taken at the surface and a prescribed atmospheric temperature TS, chosen as(2.2)where T0=15°C is a reference temperature and ΔT is the temperature difference between the southern and northern latitude of the domain. The forcing is distributed as a body forcing over the first (upper) layer of the ocean having a depth Hm.
Temperature differences in the ocean cause density differences according to(2.3)where αT is the volumetric expansion coefficient and ρ0 is a reference density. Inertia is neglected in the momentum equations owing to the small Rossby number; we use the Boussinesq and hydrostatic approximations and represent the horizontal and vertical mixing of momentum and heat by constant eddy coefficients. With r0 and Ω being the radius and angular velocity of the Earth, respectively, the governing equations for the zonal, meridional and vertical velocity u, v and w, respectively, the dynamic pressure p (the hydrostatic part has been subtracted) and the temperature T=T*−T0 become(2.4a)(2.4b)(2.4c)(2.4d)(2.4e)where is a continuous approximation of the Heaviside function; g is the gravitational acceleration; and is the surface adjustment time scale of heat, where Cp is the constant heat capacity. In these equations, AH and AV are the horizontal and vertical momentum (eddy) viscosity, and KH and KV are the horizontal and vertical (eddy) diffusivity of heat, respectively. In addition, the operators in the above equations are defined as
Slip and zero heat flux conditions are assumed at the bottom boundary, while at all lateral boundaries no-slip and zero heat flux conditions are applied. As the forcing is represented as a body force over the first layer, slip and zero heat flux conditions apply at the ocean surface. Hence, the boundary conditions are(2.5a)(2.5b)(2.5c)The parameters for the standard case are the same as in typical large-scale low-resolution ocean general circulation models and their values are listed in table 1.
3. The AMO mode
To determine whether Hopf bifurcations occur within the minimal primitive equation model of §2, we use methods from numerical bifurcation theory (Dijkstra 2005). First, the governing equations (2.4a)–(2.4e) and boundary conditions (2.5a)–(2.5c) are discretized on an Arakawa B-grid using central spatial differences. In the results of this section we use a horizontal resolution of 4°. An equidistant grid with 16 levels is used in the vertical so that the first layer thickness Hm=250 m. The discretized system of equation can be written in the form (1.3) and a 16×16×16 grid with five unknowns per point (u, v, w, p and T) leads to a dynamical system of dimension (the number of degrees of freedom) 20 480.
The steady equations of the form (1.3) are solved using a pseudo-arclength continuation method (Keller 1977). As the primary control parameter μ, we choose the equator-to-pole temperature difference ΔT. For every value of ΔT we calculate a steady solution of the minimal model under the restoring flux Qrest in (2.1). For each steady flow pattern the maximum of the meridional overturning streamfunction (ψM) is calculated and plotted against ΔT in figure 2a. The meridional overturning streamfunction for ΔT=20°C is plotted in figure 2b. The maximum of ψM occurs at approximately 55° N and the amplitude is approximately 16 Sv.
Next we diagnose the ocean–atmosphere heat flux Qpres of each of the steady solutions and compute the linear stability of the steady solution under the heat flux Qpres (where the subscript refers to ‘prescribed’). To determine the linear stability we solve for the ‘most dangerous’ modes of the problem (1.4), i.e. those with the real part closest to the imaginary axis and order the eigenvalues σ=σr+iσi according to the magnitude of their real part σr (the growth factor). The growth rate and period of the mode with the largest growth factor are plotted against ΔT in figure 3.
For ΔT=20°C the AMO mode has a positive growth factor and hence the background state, of which the meridional overturning streamfunction was shown in figure 2b, is unstable to the AMO mode. The period of the AMO mode is approximately 67 years at ΔT=20°C, and it decreases with increasing ΔT (figure 3). From figure 3, we also see that the growth factor of the AMO mode decreases strongly with decreasing ΔT and becomes negative near ΔTc≈4°C, where the Hopf bifurcation occurs. For ΔT<ΔTc the steady states are therefore linearly stable under the prescribed flux Qpres.
It was shown in Dijkstra (2006) that, for small ΔT, the angular frequency of the AMO mode becomes zero and the complex conjugate pair of eigenvalues splits up into two real eigenvalues. The paths of the two different modes can be followed to the ΔT=0 limit, where the eigensolutions connect to those of the diffusion operator of the temperature equation, called SST modes in Dijkstra (2006). These SST modes can be ordered according to their zonal, meridional and vertical wavenumber, and it was found that the AMO mode connects to the (0, 0, 1) SST mode and the (1, 0, 0) SST mode at ΔT=0.
For each eigenvalue σ associated with the AMO mode, there is a corresponding eigenvector X=Xr+iXi according to (1.4). In figure 4, the SST field of the real part of the eigenvector (Xr) of the AMO mode is plotted for ΔT=4°C (near a Hopf bifurcation) and ΔT=20°C. A comparison of the pattern in figure 4b and the one in fig. 4d of te Raa & Dijkstra (2002) demonstrates that the AMO mode here is the multidecadal mode as described in detail in te Raa & Dijkstra (2002). With increasing ΔT, the pattern becomes more localized in the northwestern part of the basin.
The physical mechanism of propagation of the AMO mode was presented in te Raa & Dijkstra (2002). This mechanism holds for every ΔT for which an oscillatory AMO mode is present (cf. figure 4). A slight generalization (compared with that in te Raa & Dijkstra (2002)) of this mechanism is provided with the help of figure 5. A warm anomaly in the north-central part of the basin causes a positive meridional perturbation temperature gradient, which induces (via the thermal wind balance) a negative zonal surface flow (figure 5a). The anomalous anticyclonic circulation around the warm anomaly causes southward (northward) advection of cold (warm) water to the east (west) of the anomaly, resulting in the westward phase propagation of the warm anomaly. Owing to this westward propagation, the zonal perturbation temperature gradient becomes negative, inducing a negative surface meridional flow (figure 5b). The resulting upwelling (downwelling) perturbations along the northern (southern) boundary cause a negative meridional perturbation temperature gradient, inducing a positive zonal surface flow, and the second half of the oscillation starts. The crucial elements in this oscillation mechanism are the phase difference between the zonal and meridional surface flow perturbations, and the westward propagation of the temperature anomalies (te Raa & Dijkstra 2002). The presence of salinity does not essentially change this mechanism; density anomalies will take over the role of temperature anomalies in the above description.
4. Effects of additive noise
To perform transient flow computations, we use v. 3.1 of the GFDL Modular Ocean Model (MOM; Pacanowski & Griffies 2000) on the same domain and with the same forcing, resolution, boundary conditions and parameters as provided in §2. The only difference with the results in §3 is that here a stretched grid with 16 layers is used in the vertical so that the first four layers have a thickness of Hm=50 m, with the thickness then increasing to 583 m in the lowest level. The internal mode in MOM has a time step of 1 day and the external mode has a time step of 225 s. Patterns of variability are analysed using the MSSA toolkit (Ghil et al. 2002).
Restoring and prescribed flux conditions are the two limits of atmospheric damping of SST anomalies, and, to study what happens between these limits, a new general boundary condition for the surface heat flux (QD) is chosen as(4.1)where Qrest is the same as in (2.1) and Qpres is the diagnosed heat flux of each steady state. A value of γ representative of the real ocean can be estimated by examining the damping time scale of SST in the upper layer ocean. Under the forcing (4.1), the damping time scale τT is defined as(4.2)Using variables from the model (ρ=1×103 kg m−3, Cp=4.2×103 J kg−1 K−1 and Hm=50 m) and λT=20 Wm−2 K−1 and τT=30 days as a representative mid-latitude value (Barsugli & Battisti 1998) gives a value of γ≈0.75.
When γ decreases from γ=1 (prescribed flux conditions), we expect that the growth factor of the AMO mode (as in figure 4a) will decrease as the atmospheric damping becomes larger. As for γ=0 (restoring conditions), the growth rate of the AMO mode is negative (te Raa & Dijkstra 2003), and there must also be the same Hopf bifurcation (as in §3) somewhere between γ=0 and 1. Under the heat flux (4.1), we therefore compute equilibrium states using time integration starting from the γ=0 steady solution. For each solution obtained, the standard deviation of the SST over the box B: [46° N–62° N]×[74° W–50° W] is plotted in figure 6a. For γ<0.85, there are no oscillations and, near γc=0.85, the system undergoes the Hopf bifurcation and the multidecadal oscillations appear for γ>γc. The amplitude of the oscillations is measured by calculating the standard deviation of the box-averaged SST rather than the peak-to-peak amplitude, so that comparisons can be made with later simulations where the variability is not regular. The oscillations have periods decreasing from 53 yr at γ=0.85 to 45 yr at γ=1, and so the change in period with γ is much smaller than that with ΔT. The first two EOFs of the SST field, which together explain 92.0% of the variance, are shown for γ=0.9 in figure 6b. We can see that the variance in temperature at the sea surface is concentrated in the northwest region of the basin similar to that of the AMO mode (figure 4b).
Next, we consider the effect of spatial and temporal coherence in the noise forcing under the conditions γ<γc by comparing the responses of the model under the following two heat fluxes:(4.3a)(4.3b)In QW, λ is the amplitude of the noise and Zij is a normally distributed random variable that takes on a different value at each grid point (i, j) in space at each time step t. The noise in QW is thus uncorrelated in both space and time. In , and are the grid variables in the x and y directions. Furthermore, Zm(t) is a normally distributed random variable, where m indicates the number of days that this variable is persistent. The spatial pattern in (4.3b) is chosen as a rough approximation to variations in atmospheric heat fluxes seen over the North Atlantic (Cayan 1992), such as those associated with the North Atlantic Oscillation. In both QW and each , the amplitude (λ) of the noise was taken to be 10% of the difference between the minimum and maximum over the basin of the prescribed heat flux Qpres, which is approximately 20 W m−2.
In figure 7, the spectra for the case of γ=0.8 with the addition of five different types of noise are shown. Although the noise added to the system has no preferred frequency, the spectrum shows a large peak at multidecadal frequencies. This is in contrast to the case of γ=0.8 in the absence of noise, where neither the temperature nor the overturning strength vary at all. It can also be clearly seen that both the spatial and the temporal correlations of the noise increase the height and breadth of the multidecadal peak. The multidecadal peak increases as the time scale of the persistence of the forcing increases. When the spatial coherence is removed (so that the noise added to each grid point is independent), the temporal coherence still causes large variations in temperature, but the power at multidecadal frequencies is greatly reduced (not shown).
When γ is decreased below the critical value γc, noise dominates the patterns seen in the MSSA. However, if the data are low-pass filtered to allow periods from 30 to 100 years, then the patterns of multidecadal variability can be seen. Figure 8a–d shows the first four EOFs of SST for the case, with the EOFs accounting for 23.8, 23.2, 8.2 and 8.0% of the variance, respectively. Since the eigenvalues of these EOFs are so closely paired, we must take into account North's rule of thumb (North et al. 1982), which shows that the smaller the difference between the eigenvalues of two EOFs, the larger the error. In this case, the model integrations were long enough (2000 yr) that the approximations for the typical error found using North's rule of thumb are small, giving some confidence in the interpretation of the EOFs. Signals from the sinusoidal spatial pattern of the noise forcing are evident, particularly in EOF 4 as well as in the southern part of the basin in EOFs 1 and 3. The EOFs, however, still display the pattern of the AMO mode that has been excited by the noise forcing.
Figure 9 shows the effect of the different noise forcing on the standard deviation of the SST in the box B. For values of γ>γc, the noise has only a small effect. By contrast, for values of γ near and below γc, the noise causes the surface temperature to increase substantially. With the flux , the largest amplitude of the variability is achieved, showing that the spatial and temporal coherence in the noise are important to set the amplitude of the multidecadal variability. In addition, the amplitude of the variability versus γ, as shown in figure 9, is qualitatively very similar to that of the stochastic Hopf bifurcation in figure 1, showing that the AMO mode is excited by the atmospheric noise. The mechanism of this excitation is outside the scope of this paper and is examined in detail in Frankcombe et al. (submitted).
5. Summary and discussion
The aim of this paper was to show that a similar dynamical systems framework can be formulated for the AMO as for ENSO. The minimal primitive equation model, as presented in §2, takes the same role in the AMO theory as the ZC model in the ENSO theory. In the ZC model the ocean–atmosphere coupling strength μ serves as the main control parameter, while in the AMO theory the atmospheric damping time scale of SST anomalies, here mimicked by the parameter γ, has that role.
In both the minimal models of the ENSO and the AMO there is a normal mode, the most dangerous mode (having the largest growth factor), which is able to destabilize the background state. The nature of the ENSO mode that destabilizes the Pacific mean state at sufficiently large coupling strength is a merger between an ocean basin mode and a stationary SST mode. The mechanism of the ENSO mode propagation and time scale is known to be related to the dominant feedback mechanism (thermocline, upwelling and zonal advection) and the equatorial wave propagation (Neelin et al. 1998). In the AMO model, the most dangerous normal mode is called the AMO mode, and it results from a merger of two SST modes at small ΔT (Dijkstra 2006). The propagation of the AMO mode is determined by a thermal wind response to a propagating temperature (or density) anomaly and the multidecadal time scale is set by the east–west propagation time of the temperature (or density) anomalies. The oscillation can be characterized by an out-of-phase response of the meridional and zonal overturning anomalies, as shown in figure 5.
It is not known whether the tropical Pacific climate state is near a Hopf bifurcation (Fedorov & Philander 2000). Owing to slow variations of this background state, it is possible that one ENSO event could occur in the supercritical regime (with less impact of noise) and the next in the subcritical regime (with noise controlling its amplitude). Similarly, it is not known whether the atmospheric damping time scale of temperature anomalies in the North Atlantic (in this paper mimicked by γ) induces a positive or negative growth factor of the AMO mode in the deterministic case, as this depends also on the background state.
When the additive noise is added, the ENSO mode can be excited below critical conditions (Roulston & Neelin 2000) and the ENSO variability results from a stochastic Hopf bifurcation as introduced in §1. Also for the AMO, the results here make a case that a stochastic Hopf bifurcation occurs. Oscillatory variability with an amplitude depending on the noise arises for values of γ≤γc, while for γ>γc the variability does not differ much from the deterministic case. The presence of the noise forcing continuously excites the variability, and a spectrum shows that the variability has the same multidecadal period as in the cases where γ>γc. Both the spatial and temporal coherence in the random part of the heat flux forcing are important to excite the multidecadal variability to a reasonable amplitude.
When the minimal model is run at a small value of γ such that the background state is very stable, each of the noise forcing is only able to cause very small variability (not shown). Hence the presence of the AMO mode and the occurrence of the Hopf bifurcation certainly play a central role in the amplitude of the multidecadal variability. Other possible mechanisms such as a passive response of the ocean on the atmospheric noise (Hasselmann 1976) or the extension by Saravanan & McWilliams (1998), where the effect of horizontal advection leads to a preferred time scale, are therefore less probable. It is also interesting that, in the noise-forced cases, one sees normal mode patterns in the variability instead of non-normal mode patterns (Farrell & Ioannou 1996). It is probable that the time scale on which non-normal modes grow is much faster than the typical time scale of the variability. Without knowing the non-normal modes for the minimal primitive equations model, it is difficult to assess their role in the multidecadal variability.
These results are for the minimal model, but how about more realistic models, i.e. extensions of the minimal model? Until now, only the AMO mode and the periodic oscillations under prescribed flux conditions have been studied in extensions of the minimal model. While the continental shape is irrelevant for the ENSO, the shape of the continents is essential for the deformation of the AMO mode into a pattern that resembles the patterns obtained in coupled climate models (Delworth & Greatbatch 2000) and observations (Dijkstra et al. 2006). Results in idealized models with two basins showed that the AMO mode is localized in the sinking regions of the global thermohaline flow (von der Heydt & Dijkstra 2007). This indicates that the AMO mode is unique to the North Atlantic.
In summary, the results provided here indicate that a basic dynamical systems framework of the AMO can be formulated and that this framework is similar to that in the ENSO theory. The central element is that the excitation of the AMO mode by the atmospheric noise is causing the variability associated with the AMO. Although it is not easy to falsify this theory by the instrumental record, as the time scale is rather long in relation to the available data, we hope that the basic ideas will stimulate further analysis of model results and observations.
This work was funded by the Dutch Science Foundation (Earth and Life Sciences) through project ALW854.00.037 (L.M.F. and H.A.D.) and a VENI-grant (A.S.vdH.).
One contribution of 12 to a Theme Issue ‘Stochastic physics and climate modelling’.
- © 2008 The Royal Society