# Low-frequency variability and heat transport in a low-order nonlinear coupled ocean-atmosphere model

## Abstract

We formulate and study a low-order nonlinear coupled ocean–atmosphere model with an emphasis on the impact of radiative and heat fluxes and of the frictional coupling between the two components. This model version extends a previous 24-variable version by adding a dynamical equation for the passive advection of temperature in the ocean, together with an energy balance model.

The bifurcation analysis and the numerical integration of the model reveal the presence of low-frequency variability (LFV) concentrated on and near a long-periodic, attracting orbit. This orbit combines atmospheric and oceanic modes, and it arises for large values of the meridional gradient of radiative input and of frictional coupling. Chaotic behavior develops around this orbit as it loses its stability; this behavior is still dominated by the LFV on decadal and multi-decadal time scales that is typical of oceanic processes. Atmospheric diagnostics also reveals the presence of predominant low- and high-pressure zones, as well as of a subtropical jet; these features recall realistic climatological properties of the oceanic atmosphere.

Finally, a predictability analysis is performed. Once the decadal-scale periodic orbits develop, the coupled system’s short-term instabilities — as measured by its Lyapunov exponents — are drastically reduced, indicating the ocean’s stabilizing role on the atmospheric dynamics. On decadal time scales, the recurrence of the solution in a certain region of the invariant subspace associated with slow modes displays some extended predictability, as reflected by the oscillatory behavior of the error for the atmospheric variables at long lead times.

###### keywords:

Extended-range predictability, Low-frequency variability (LFV), Low-order modeling, Lyapunov instability, Ocean-atmosphere coupling, Slow periodic orbit.###### Msc:

[2014] 00-01, 99-00^{1}

[mycorrespondingauthor]Stéphane Vannitsem

## 1 Introduction and motivation

The variability at annual, interannual and decadal time scales of the coupled ocean–atmosphere system is currently a central concern in improving extended-range weather and climate forecasts. The oceans’ long-term variability and their interaction with the atmosphere has been extensively explored in the Tropical Pacific, due to the climatological importance of the El Niño–Southern Oscillation (ENSO) phenomenon (e.g. (14); (51); (64)). The ocean and the atmosphere also interact in the mid-latitudes through both wind stress and buoyancy fluxes (13); (66)), although the impact of this interaction on the long-term variability of the atmosphere is still controversial ((4); (35)).

On physical grounds, interactions between the two components of the coupled climate system in mid-latitudes are obviously essential to its functioning on multiple time scales. The main direction of the coupling, however, is a matter of debate: Is the slower ocean responding to the wind stress forcing in an essentially passive way (19) — i.e., is its feedback to the atmosphere too weak to qualitatively modify the dynamics of a stand-alone atmosphere — while, at the same time, playing the role of a heat bath that provides boundary forcing for the atmosphere (55)? Or is the ocean playing a more active role in atmospheric dynamics Feliks et al. (2004); Feliks et al. (2007); (28)? Marshall et al. (2001), for instance, discuss these questions in detail.

From a dynamical point of view, one possible answer to these questions translates into a search for the presence of coupled modes between the ocean and the atmosphere, such as found for instance in observational data (9). In the context of coupled global climate models, the answer, however, differs from one investigation to another or from a modeling setup to another; this answer depends, to a large extent, on whether the forcing is generated by deterministic or stochastic, linear or nonlinear processes Delworth and Mann (2000); Kravtsov et al. (2008); (39).

Part of our motivation is that a low-frequency nonlinear oscillation has been found in intermediate-order, coupled nonlinear ocean-atmosphere models. This coupled mode operates by triggering the displacement of the atmospheric jet position Kravtsov et al. (2007); (67). The physical origin of this coupled mode is, however, not fully understood as yet, nor is its existence generally agreed upon.

To understand the qualitative behavior of the coupled ocean–atmosphere dynamics, several low-order models have been developed. Such models allow one to isolate the essential processes believed to be at play in a specific problem at hand, by using as building blocks only the minimal ingredients describing the dynamics.

This approach — originating in the works of Saltzman (57) and Lorenz (40) and, from then on, in the development and applications of nonlinear dynamics to the environmental sciences — attempts to reduce complicated dynamics to its essential features. It has been quite successful, so far, in increasing our understanding of the dynamics of the ocean alone, in particular the multimodality of the thermohaline circulation (13); (63), as well as the development and variability of the mid-latitude oceanic gyres (28); Pierini (2012); (59); Veronis (1963).

To the best of our knowledge, it is Lorenz (42) who applied this approach first to the coupled system and developed a pseudo-spectral, low-order model based on the primitive equations for the atmosphere, coupled to an ocean heat bath. Vannitsem (72) used this model to evaluate the impact of climate changes on model error biases. Nese and Dutton (47) extended the model by adding an ocean dynamics similar to the one proposed by Veronis (1963) and based on four dominant ocean modes. These authors found that accounting for the heat transport helps increase the coupled model’s predictability.

Other minimal-order coupled models (56); Van Veen (2003) allowed for the possibility of the ocean’s developing a thermohaline circulation. In the latter work, Van Veen (2003) performed a bifurcation analysis and showed the active role of the ocean in setting up the dynamics when close to the bifurcation points of the atmospheric model. While Deremble et al. (2012) did not consider full coupling, they did point out some of the similarities between the bifurcation trees of the atmospheric and oceanic dynamics in a mid-latitude setting.

More recently, Vannitsem and colleagues (73); (74) have developed two coupled model versions based on a quasi-geostrophic atmospheric model proposed by Charney and Straus (1980) and further studied in (53); (54), and on the ocean model of Pierini (2012). The coupling in both of these versions was based solely on a mechanical transfer of momentum via wind friction. In particular, these authors showed that their coupled model displays decadal variability within the ocean, similar to that found in idealized, intermediate (28); (62) and low-order models (5); (46); (59). Furthermore, the investigation of this coupled model’s stability properties showed that the momentum transfer coupling between the two components contributes to an increase of the flow’s instability, in terms of both the magnitude of the positive Lyapunov exponents and their number.

This model (73); (74) is, however, missing an important ingredient of the coupled dynamics, namely the energy balance between the ocean and the atmosphere. The present work proposes a new model version, in which the thermal forcing affecting only the atmosphere is replaced by an energy balance scheme (3); Chen and Ghil (1996); Deremble et al. (2012). It will allow us to disentangle the respective roles of the heat and radiative fluxes through the ocean surface vs. the transport of heat within the ocean in affecting the coupled system and its predictability.

The model equations are described in Section 2 , for the dynamics (73); (74) and thermodynamics, respectively; they are reduced to a low-order system in Section 2.5. In Section 3, a bifurcation analysis of the basic solutions is first performed (Section 3.1); it reveals the presence of low-frequency variability (LFV) in the form of a set of long-periodic, attracting orbits that couple the dynamical modes of the ocean and the atmosphere. The model dynamics is further explored through the analysis of the climatological properties of the solutions in Section 3.2, while the dependence of the decadal-scale orbits on model parameters and their predictability are studied in Sections 3.3 and 3.4, respectively. A summary of the results and future prospects are then provided in Section 4.

## 2 The model equations

### 2.1 Equations of motion for the atmosphere

The atmospheric model is based on the vorticity equations of a two-layer, quasi-geostrophic flow defined on a -plane (50); (66). The equations in pressure coordinates are

(1) | |||||

here and are the streamfunction fields at 250 and 750 hPa, respectively, and is the vertical velocity, is the Coriolis parameter at latitude , with its meridional gradient there.

The coefficients and multiply the surface friction term and the internal friction between the layers, respectively, while hPa is the pressure difference between the two atmospheric layers. An additional term has been introduced in this system in order to account for the presence of a surface boundary velocity of the oceanic flow defined by (see the next subsection). This would correspond to the Ekman pumping on a moving surface and is the mechanical contribution of the interaction between the ocean and the atmosphere.

### 2.2 Equations of motion for the ocean

The ocean model is based on the reduced-gravity, quasi-geostrophic shallow-water model on a -plane (e.g., (26); Pierini (2012); (66)):

(2) |

Here is the streamfunction in the model ocean’s upper, active layer, which has density , depth , and lies over a quiescent deep layer with density ; is referred to as the reduced gravity felt by the fluid in the active layer, and is the reduced Rossby deformation radius. The friction coefficient at the bottom of the active layer is , and is the vertical component of the curl of the wind stress.

Usually, in low-order ocean modeling, this curl is prescribed as an idealized profile, zero along a latitude placed at or near 45N, and antisymmetric about this latitude (28); (59). In the present work, this forcing is provided by the wind stress generated by the atmospheric component of the coupled system.

Assuming that the wind stress is given by — where are the horizontal components of the geostrophic wind, and the corresponding components of the geostrophic currents in the ocean — one gets

(3) |

Here the wind stress is proportional to the relative velocity between the flow in the ocean’s upper layer and the wind in the lower atmospheric layer. The drag coefficient characterizes the strength of the mechanical coupling between the ocean and the atmosphere and will be a key bifurcation parameter in the present work.

### 2.3 Ocean temperature equation

We assume here that temperature is a passive scalar transported by the ocean currents, but the oceanic temperature field displays strong interactions with the atmospheric temperature through radiative and heat exchanges (3); Deremble et al. (2012). Under these assumptions, the evolution equation for the ocean temperature is

(4) |

with

(5) |

In Eqs. (4) and (5) above, is the net radiative flux at the ocean surface. This flux contains three terms: the shortwave radiation entering the ocean; the outgoing longwave radiation flux ; and the longwave radiation flux re-emitted to the ocean; here is the emissivity of the atmosphere and is the Stefan-Boltzmann constant. The parameter is the heat capacity of the ocean, and is the heat transfer coefficient between the ocean and the atmosphere that combines both the latent and sensible heat fluxes. We assume that this combined heat transfer is proportional to the temperature difference between the atmosphere and the ocean. The typical values of the above parameters are provided in Table 1.

### 2.4 Atmospheric temperature equation

Let us now consider the equation for the temperature of the baroclinic atmosphere presented in Section 2.1. As in the ocean, a radiative and heat flux scheme is incorporated reflecting the exchanges of energy between the ocean, the atmosphere and outer space (3):

(6) |

with

(7) |

In Eq. 6, is the atmospheric barotropic streamfunction, is the net radiative flux in the atmosphere, the gas constant, the vertical velocity in pressure coordinates, and

is the static stability, with the pressure, the air density, and the specific heat at constant pressure; here is taken to be constant. Note also that, thanks to the hydrostatic relation in pressure coordinates and to the ideal gas relation , the atmospheric temperature can be expressed as .

As for the radiative flux in the ocean, in Eqs. (4) and (5), the net radiative flux within the atmosphere is composed of three terms: the ingoing flux of radiative energy effectively absorbed; the outgoing flux re-emitted to the ocean and to space; and the shortwave radiative flux absorbed directly by the atmosphere. Note that, in writing Eqs. (5) and (7), we assume that the rates of radiative emission and absorption are equal, i.e that emissivity is equal to absorptivity. This assumption is strictly valid only at thermodynamic equilibrium, but it can be safely applied to systems in local thermodynamic equilibrium, like the lower atmosphere, in which molecular collisional processes dominate the radiative processes (27).

### 2.5 Low-order model formulation

In order to build a low-order model version, we follow the truncated Fourier expansion approach in Charney and Straus (1980); (28); (40); Pierini (2012); (53); (57); (59); (74); Veronis (1963), i.e., the model fields are developed in a series of basis functions and truncated at a minimal number of modes that still captures key features of the observed behavior. Both linear and nonlinear terms in the equations of motion are then projected onto the phase subspace spanned by the modes retained, by using an appropriate scalar product.

In the thermodynamic equations introduced in Sections 2.3 and 2.4, quartic terms appear in the radiative fluxes. To overcome this problem, we will take advantage of the small amplitude of temperature anomalies, as compared with the global mean, in order to linearize these terms. The details of this linearization are described in B.

For the model’s closed ocean basin, we use only sine functions, in order to avoid fluxes through the boundaries. For the atmosphere, no-flux boundaries are assumed in the meridional direction and periodicity in the longitudinal direction. Hence both sine and cosine functions are allowed in the longitudinal direction, while we use cosine modes only in the meridional direction, as discussed for instance in (53). The radiative fluxes are taken as proportional to a cosine function of latitude.

To obtain the coupled-model equations for the atmosphere, we combine the truncated temperature equations (23) of C with the equations for the atmospheric streamfunction perturbation that result from the projection of the dynamical equations (1) onto the modes in the set (C). As in (53), one obtains a set of ordinary differential equations (ODEs) for the atmosphere.

Similarly, combining the projection of the streamfunction field for the ocean onto the modes of (C) with the temperature field of (22) leads to a set of ODEs, cf. C. The full coupled model is thus based on a set of ODEs. The time-dependent solutions of this system are obtained by numerical integration, using a second-order Heun method with a fixed time step time units. Several higher-accuracy methods have been tested but without affecting substantially the results reported herein.

In our analysis of the dynamics of the coupled system, we will mostly focus on three parameters, given in dimensional units, in order to ascertain the typical values that give rise to remarkable behavior. These three parameters are [Wm], [s and [WmK]: they correspond to the meridional variation in the radiative input from the Sun; the strength of the coupling between the ocean and the atmosphere, as an inverse of a response time scale; and the intensity of the heat fluxes, respectively.

## 3 Dynamics of the coupled model

### 3.1 Bifurcation diagram

We start our analysis of the coupled model’s dynamics by constructing its bifurcation diagram, cf. (13); (21); (60) and references therein. The AUTO software (12) is used to follow solution branches by pseudo-arclength continuation and detect local bifurcations. The bifurcations will be explored in the two-dimensional parameter space of and , with and WmK fixed.

In constructing the model’s bifurcation diagram, we first fix the mechanical friction coefficient between the atmosphere and the ocean at a value of s, and vary the parameter that scales the energy absorption by the ocean between 0 and 400. The norm plotted in Fig. 1a summarizes the successive bifurcations of the solutions. We define the -norm of a solution by

for a fixed point , and | |

for a periodic orbit with , |

where is the period.

When the system has only one fixed point, which remains stable as increases, until a first Hopf bifurcation occurs at Wm (Fig. 1a); it is clearly subcritical.

The periodic orbits along this branch have very long periods, of roughly 21 years. This “slow” branch encounters a fold at Wm and, from there on, it stabilizes and extends forward as the parameter increases further. Such a fold is often associated with a “global” Hopf bifurcation, in the terminology of (23); see also ((21), Ch. 11) and references therein. The stable periodic branch loses its stability again at Wm, according to the values of the Lyapunov exponents computed in Section 3.4.

Interestingly, the long-periodic branch is no longer present in Fig. 1b, i.e. at . This essential difference between panels (a) and (b) of Fig. 1 indicates that the frictional coupling between ocean and atmosphere is at the origin of the development of the long-periodic oscillation. In other words, wind friction is, in the present coupled model, the main physical process triggering the LFV development within both the atmosphere and the ocean.

After the first Hopf bifurcation, the fixed-point branch undergoes another Hopf bifurcation, at Wm. The periodic orbits that emanate from this bifurcation have short periods, of 10–15 days, and are unstable. According to the AUTO software, this “fast” branch possesses many torus bifurcations, also referred to as secondary Hopf bifurcations. Indeed, depending on the value of , various stable quasi-periodic solutions are found in the neighborhood of this branch: these solutions are characterized by the coexistence of one long period, of about 22 years, and a short one, of about 5 to 20 days. Concomitantly, these solutions are associated with a weak transport, whose intensity is much smaller than the transport associated with the slow branch.

As shown in Fig. 1b, the fast branch is already present in the bifurcation diagram when the ocean and the atmosphere are not coupled mechanically, i.e. at . For this value of , the fast branch is stable until a torus bifurcation occurs at Wm. The quasi-periodic solutions found after this bifurcation do not possess the above-mentioned long periodicity, only the short ones. These findings indicate the purely atmospheric origin of the latter, short-periodic oscillations. Indeed, the oceanic transport vanishes for the quasi-periodic solutions near this secondary Hopf bifurcation, . Interestingly, along the fast branch and for the tori found nearby, all the atmospheric modes and a few of the ocean temperature modes — namely , with — are oscillating on the fast time scale, while the other modes remain stationary.

Oscillatory atmospheric modes having periodicities that are longer than the typical life time of extratropical storms but shorter than a season, i.e. roughly of days, are called intraseasonal and have been described and discussed at some length in the dynamic-meteorology literature, especially in the Northern Hemisphere extratropics (21); (22). Their impact on similarly fast oceanic modes in this coupled model is interesting but not entirely surprising.

The loci of the two Hopf bifurcations that lead to the slow and the fast periodicities intersect in a Hopf–Hopf bifurcation, as seen in the regime diagram of Fig. 2. The locus associated with the slow branch finally terminates in a Bogdanov-Takens bifurcation. A locus of torus bifurcations emanates from the Hopf–Hopf bifurcation and accounts for the bifurcations of this type found along the fast branch in Fig. 1a. Another locus of torus bifurcations accounts for the bifurcation that is already present when ; see discussion of Fig. 1b above. The locus of the fold of the long-periodic orbit is depicted in red in Fig. 2 and we see that it merges with the curve of the Hopf bifurcation of the long-periodic orbit in a generalized Hopf bifurcation at Wm.

This bifurcation diagram provides only a partial view of our coupled model’s rich dynamics, due to two factors: first, several additional parameters may trigger further qualitative changes in system behavior, like the heat flux parameter ; and second, other bifurcation types, like global bifurcations or blue-sky catastrophes, might be present ((36); (60)). The rich behavior presented in Sections 3.2 and 3.4 below may arise, at least in part, from the presence of such nonlocal bifurcations.

Nevertheless, the local bifurcation diagram already reflects the wide variety of solutions generated by this relatively simple model, and in particular the presence of a long-periodic solution associated with the coupling of the atmosphere with the ocean. We must emphasize that this bidecadal periodicity is not only a property of the ocean but also of the dominant dynamical modes of the atmosphere, as illustrated in Fig. 3. In this figure, we plot a three-dimensional (3-D) projection of the long-periodic orbits onto the subspace spanned by the modes , for and for several -values.

All the periodic orbits plotted in Fig. 3 involve not just the three variables shown, but a total of variables — namely for , for , for , and for — while the other variables are equal to 0. If the latter variables are set to 0 initially, they remain equal to 0 as the flow evolves. Hence the subspace of the 17 variables listed above is invariant under the phase space flow induced by the 36 ODEs that govern our coupled model. Moreover, all the orbits we computed within this subspace were dominated by slow motions; hence the use of the subscript ‘s’ for “slow” and of ‘f’ for “fast” with respect to the 19 other variables.

The orbit shown for Wm (blue curve in Fig. 3) is not only periodic but also attracting, i.e., it is a genuine, stable limit cycle. The other periodic orbits plotted in the figure, across the range of parameters Wm, while unstable for arbitrary perturbations, are actually stable to perturbations within the slow subspace, even for very large values of .

Perturbations introduced in the complementary 19-dimensional subspace — for values of larger than about Wm — do not die out, but saturate fairly quickly ; see Section 3.4, especially Fig. 15 there. It is reasonable, therefore, to conjecture that the limit cycles in Fig. 3 form the “backbone” of a strange attractor. The dimension of this attractor will be examined in Sections 3.3 and 3.4 below, while the precise meaning we attach to the term backbone will be discussed in Section 4.

We will see in fact later that the LFV associated with the long-periodic orbits in Fig. 3 — whether stable or not — is still present when chaotic solutions develop beyond a certain range of parameter values (). The decadal-scale limit cycles thus continue to organize the coupled-system dynamics, well beyond this range, and the chaotic solutions still live close to them.

### 3.2 Climatological properties for large values of

Before discussing in further detail the model dynamics around the long-periodic orbits shown in Fig. (3), let us now focus briefly on the climatological properties of typical solutions generated by the system. The numerical code is the same second order code already mentioned in the previous section.

Figure 4 displays the mean fields for the ocean and the atmosphere for Wm and s. The ocean is displaying in panel (c) the double-gyre dynamics that is well known from ocean models with prescribed, time-independent wind stress ((13); (28); (60)), but now in a genuinely coupled ocean–atmosphere model and thus including the temperature field in panel (d). This double gyre is transporting — in the present, coupled model — heat toward the pole, and it thus reduces the heat contrast affecting the atmosphere through its interaction with the ocean, as is the case in the observations ((13)) and in much more detailed models (L’Hévéder et al., 2014, and references therein).

In addition, a clear high- and low-pressure dipole is appearing in the atmospheric streamfunction and temperature fields of Figs. 4(a,b), with a low-pressure center in higher latitudes, north of 45N, and a high-pressure center further south. The atmospheric dipole is present whatever the value of , i.e., it is independent of the nature of the dynamics, whether stationary, periodic or chaotic. This feature is also present when , indicating that it is induced by the model’s radiative scheme (not shown).

For all parameter values examined, a mean jet is forming between the two pressure centers. In the present, coupled model, the two pressure centers, and hence the jet between the two, are strengthened by the localized heat transfer between the atmosphere and the ocean.

The mean dipole and the jet between the two “semi-permanent centers of action” are reminiscent of similar features in the real atmosphere over the North Atlantic and North Pacific (Lorenz, 1967, and references therein). For instance, the Icelandic low and the Azores high control the position and intensity of the mid-latitude jet affecting the Atlantic and Western Europe.

The intensity of the low- and high-pressure centers also depends strongly in our coupled model on the degree of mechanical coupling between the ocean and the atmosphere, as shown in Fig. 5. Furthermore, the increase of the coupling reduces the temperature contrast within the ocean, which in turn reduces the one within the atmosphere; compare panels (b) and (d) of Figs. 4 and 5. We suspect, therefore, that the predictability of the system will be highly affected by the ocean transport, since a drastic reduction of the meridional temperature gradient is experienced. We investigate this aspect in Section 3.4, in which Lyapunov exponents are computed.

Feliks et al. (2004, 2007) have shown that the mid-latitude jet is much stronger yet in the presence of a sharp sea surface temperature (SST) front, like the Gulf Stream or the Kuroshio. Such a sharp front is not possible in our intermediate, low-resolution coupled model, nor in typical general circulation models (GCMs) used until recently in climate change studies. When such a front is introduced into a GCM, the stronger jets predicted in (Feliks et al., 2004, 2007) do appear ((4); Minobe et al. (2008)).

### 3.3 Decadal-scale dynamics and its dependence on coupling

In order to better understand the changes in model dynamics that occur when the mechanical coupling is increased, we have represented in Fig. 6 the projection of the trajectories onto the 3-D subspace spanned by the modes , as in Fig. 3. For small values of , the solutions of the coupled model are well localized around small values of ; see red and green trajectories in panels (a) and (b). When is increased, the ocean dynamics becomes more vigorous, the transport of heat toward the pole increases, and the amplitude of the LFV in this subspace increases. This increase induces the high-amplitude cycles displayed in Figs. 6a and b.

Along these large limit cycles, when the meridional temperature gradients are intense, i.e. when is large, the transport is strong, in order to reduce the gradient. At the same time, the heat transfer out of the ocean induces a high variability within the atmosphere. Once the amplitude of the temperature gradient associated with decreases, the oceanic transport becomes less intense and the atmospheric dynamics more quiescent.

Another important and recently much-debated question (e.g., Delworth and Mann (2000)) is the role played by the heat fluxes in the coupling between the ocean and the atmosphere. The effect of on the coupled model’s LFV is shown in Fig. 7, which shows the same type of trajectories as in Fig. 6, but at different - rather than -values. The dynamics is drastically modified when increasing , with very chaotic trajectories at no heat flux (red) and a purely periodic trajectory for the largest value of in the figure. This solution behavior suggests that the larger the heat flux the stabler the dynamics. This aspect of our results will be further explored in the next subsection, by computing the coupled model’s Lyapunov exponents.

We saw that our coupled model successfully simulates semi-permanent highs and lows that have some similarity with the Azores High and the Icelandic Low, as well as multi-annual LFV. One might thus wonder whether it might generate spatial features similar to the North Atlantic Oscillation (NAO) and display its decadal variability, as proposed by Feliks et al. (2011) or by Kravtsov et al. (2007, 2008).

Figure 8 shows the geopotential height difference between two points of the spatial domain — located at ( and , in the model’s nondimensional coordinates — within the atmosphere, for different values of the radiative meridional gradient and of the coupling parameter . These points correspond roughly to high- and low-pressure centers in the model atmosphere.

When the solutions stay close to the previously discussed slow, long-periodic orbits, as they seem to do in Fig. 8a, a low-frequency oscillation is clearly present for which the delay between the dominant maxima in all three time series is roughly 20 years, with secondary maxima also present in-between. This succession of high-amplitude and low-amplitude maxima is associated with a slight asymmetry in the system’s attractor between negative and positive values of , an asymmetry that strongly affects the oceanic transport.

The strong and fairly smooth LFV in Fig. 8a contrasts with the considerably smaller-amplitude and much more irregular behavior apparent in Fig. 8b. The latter is essentially dominated by the atmospheric variability when the meridional heating gradient is large and the frictional coupling is small, as is the case in panel (b). This difference further illustrates the importance of the underlying slow orbits found in the coupled system, around which the dynamics organizes for a certain range of parameters.

To disentangle the LFV present in these model simulations, we applied singular-spectrum analysis (SSA; (25); (75)) to the time series in Fig. 8. For the model run with Wm and s, the SSA spectrum plotted in Fig. 9a indicates one periodicity, with a period of 11 years, that is significant at the 95% level; this periodicity reflects the succession of consecutive maxima within all three runs in Fig. 8a. A longer period, though, of roughly 22 years — which corresponds to the distance between the dominant maxima in the time series — also agrees with the complete length of the limit cycles in Figs. 3 and 7; this bidecadal period does not appear to be significant in the analysis so far.

In order to isolate this longer period, the reconstruction of the 11-year signal (24); (25) was removed and the SSA algorithm applied again on the residual time series, as shown in Fig. 9b. A 21.9-year period is now significant at the 95% level, as expected from the visual inspection of the time series. In addition, the third, fourth and fifth harmonic of this 22-year cycle — at 7.3, 5.5 and 4.4 years — are all significant at the same level, indicating the highly anharmonic nature of the bidecadal cycle, as apparent in Fig. 8a, too.

The discussion of these results is postponed to the paper’s last section.

### 3.4 Stability and predictability of the coupled model

#### Lyapunov exponents

In dynamical systems theory, the stability of a nonlinear system, like our coupled model, is usually quantified by evaluating its Lyapunov exponents. These exponents characterize the growth or decay of small differences in the system’s initial data, and thus help generalize stability analysis from linear to nonlinear differential equations; they are evaluated in the so-called tangent space of a model trajectory (37); (49).

Lyapunov exponents can be computed using orthogonal vectors in the tangent space of a single model orbit,

(8) |

here is the Lyapunov vector corresponding to the Lyapunov exponent, in decreasing algebraic-value order (29), is the dimension of the system, and is the time step used in this computation. In the present section, we compute our coupled model’s Lyapunov exponents and use them to explore the predictability properties of the flow as a function of the parameter values. To integrate the tangent linear system and rescale the Lyapunov vectors, we used the same second-order Heun scheme, with a time step of time units, as for integrating the full nonlinear model in Sections 3.2 and 3.3.

Figures 10(a, b) display the effect of the radiative-forcing parameter on the three leading Lyapunov exponents, for and , respectively. As seen in panel (a), model behavior experiences a transition from stationary — for and Wm — through periodic, from up to about Wm, and on to chaos, with two exponents becoming positive at almost the same value of .

For a larger shortwave radiative forcing within the atmosphere, (panel b), the transition occurs slightly earlier, at rather than Wm, and the amplitude of both positive exponents is somewhat larger. But the overall features of the transition are similar in the two panels. We will therefore focus in the following on the case for which .

As mentioned in Section 3.1, we evaluated the stability of the slow branch — green solid curve in Fig. 1a for stable solutions and green dashed for unstable solutions — by computing the Lyapunov exponents. To do so, we followed the algorithm outlined above, starting from initial points close to the long-periodic orbit, for each set of parameter values explored. Stable periodic solutions, i.e. genuine limit cycles, correspond to the first exponent in Fig. 10a being zero, and they lie precisely along the solid green line displayed in Fig. 1a. The detailed mechanism of destabilization of this slow limit cycle is not clear yet and it is probably associated with a global bifurcation, whose analysis is beyond the scope of the present study.

Figure 11 shows the complete Lyapunov spectra for Wm, in red for s and in blue for s. The spectra indicate that the solutions are chaotic for both values of . But clearly a drastic change occurs when one increases the mechanical coupling parameter , namely a substantial drop in the positive Lyapunov exponents. This drop is consistent with the emergent effect of slow limit cycles on their vicinity in phase space, as already discussed in connection with Figs. 1c, 6 and 7 in the previous three subsections.

The key role of these slow solutions in the coupled model’s dynamics is further illustrated in Fig. 12, in which we plotted the phase space distribution of the local Kolmogorov-Sinai entropy , given by the sum of the positive local Lyapunov exponents ,

(9) |

here is the state vector in the model’s 36-dimensional phase space, and is the number of positive local Lyapunov values. Kolmogorov-Sinai entropy, also called metric entropy, is discussed in ((61)), and local Lyapunov exponents in ((1)). The local metric entropy gives a measure of the overall divergence of trajectories in a neighborhood of a point , although local Lyapunov exponents, unlike the global ones, are not invariant under a nonlinear change of variables.

Figure 12 clearly indicates that the neighborhood of the slow periodic solutions is quite stable at the parameter values used in the figure. When the ocean transport is more intense and the meridional temperature gradient is larger than in the run used in Fig. 12, the flow is much more unstable, as the atmospheric dynamics is much more active (not shown).

Vannitsem and Nicolis (1997) have investigated the properties of local Lyapunov exponents and of the corresponding local Kolmogorov-Sinai entropy in a three-level quasi-geostrophic model with a fairly realistic climatology and variability for the extratropical atmosphere (Kravtsov et al., 2009; Marshall and Molteni, 1993). In that model — often referred to as the QG3 atmospheric model Kondrashov et al. (2004) — the values of the local exponents were highly dependent on the large-scale weather regimes ((70); (71)). In the present analysis, a similarly strong dependence on large-scale flow features is also evident, as indicated by the large variability of the Kolmogorov-Sinai entropy when moving along the slow solutions in Fig. 12.

Figure 13 illustrates the change in stability as a function of the coupling parameter . For Wm, the increase in first leads to a sharp drop in the leading Lyapunov exponents. Increasing further induces a transition to the periodic solution already illustrated in Section 3.3. For Wm the picture is slightly different, since there is no transition to a periodic solution for the set of values explored. The solution of the system is still chaotic, but it is living close to the slow periodic solutions that do exist when is large enough.

Changes in the heat flux parameter can also affect strongly the stability of coupled-model solutions. In Fig. 14, we plot the leading Lyapunov exponent as a function of for different combinations of the other two parameters, and . Clearly decreases as a function of , except for Wm and s. The decrease is sharpest for Wm and s (dotted blue curve) and weakest for Wm and s (dotted green curve). Overall, these results suggest that the heat fluxes between the ocean and the atmosphere have a stabilizing effect on the coupled model.

This stabilizing effect of the fluxes is also seen in examining the dimension of the attractor associated with the slow behavior in our coupled model. We consider the Kaplan-Yorke dimension (30), often called the Lyapunov dimension; it is given by a simple functional of the Lyapunov exponents,

where is the largest such that . and, under fairly general circumstances, it equals the information dimension. Hence, is preferable to several other ways of estimating the dimensionality of attractors ((65), Ch. 2).

The results are very instructive, indeed, for the three orbits plotted in Fig. 7 and for two additional ones, also computed at Wm and s. For the slow limit cycle at WmK (blue curve in the figure), one obviously has As decreases, one gets: for (not in the figure) , for (green curve) , for (not in the figure) , and for (i.e., no heat flux at all; red curve in the figure) . Visually, the green and red objects in Fig. 7 look successively stranger, while our Kaplan-Yorke dimension calculations indicate that the corresponding attractors, at Wm and s, have for all but the two lowest -values, namely 1.0 and 0.

#### Predictability

Finally, the long-term predictability of the atmospheric component is explored when the slow variability is starting to dominate coupled-model dynamics. The experiment consists in evaluating the growth of initial errors in forecasts that are issued for Wm and s.

The errors are introduced into all the variables of the coupled model, but we are most interested in error evolution for the atmospheric modes: their root-mean-square growth is tracked by restricting the -norm defined in Section 3.1 to these 20 modes according to

here and are the projection of the control and perturbed trajectories, respectively, onto the atmospheric variables.

An ensemble of 1000 realizations was obtained by using random initial errors sampled from a Gaussian distribution with mean 0 and variance , in nondimensional units, around two specific locations on the coupled model’s attractor. These two locations are indicated in Fig. 12 by the black filled circle and the green triangle, respectively, and the error growth curves are displayed in Fig. 15.

For short times and for the first initial state (black curve) the error grows very rapidly up to a time of the order of a few days. This contrasts with the short-time behavior of the error for the second initial state (green curve), which stays quite small and thus indicates potential for an accurate extended forecast beyond several years. The difference between the short-time error growth in the two cases reflects, of course, the difference in local stability along the attractor, as revealed by the Lyapunov exponents (not shown) and the Kolmogorov-Sinai entropy; see Fig. 12.

For long lead times, the error evolution exhibits a large-amplitude cycle, with large errors when the system is in the high-atmospheric-variability region along the attractor and small errors when the system revisits the region of low atmospheric variability; see again Fig. 12. As discussed in Section 3.3, these two regions coincide with the temperature gradient within the ocean.

The large up- and downswings in error evolution for long lead times in our coupled model stand in sharp contrast with the error evolution in certain stand-alone atmospheric models, like Lorenz’s low-order model of a moist general circulation ((42); (48)) or the already mentioned QG3 model (Marshall and Molteni, 1993; Vannitsem and Nicolis, 1997). In these atmospheric models, the error growth is predominantly monotonic — like in typical high-resolution numerical weather prediction models (29) — and saturates quickly at an asymptotic mean value after a typical time scale of the order of the inverse of the first Lyapunov exponent.

The oscillatory nature of the error evolution in the coupled ocean–atmosphere model at hand is the by-product of the modulation of the atmospheric variability by the oceanic modes, in the model regime in which the destabilized slow limit cycle still affects substantially the chaotic variability that replaces the purely periodic flow in the model’s phase space. This nearly periodic modulation of the error evolution — which is clearly visible in both the black and the green curves in Fig. 15 — implies, in turn, that the typical time scale of convergence of the probability density of the coupled model’s variables toward their asymptotic probability density must be much larger than the time scale associated with the dominant Lyapunov exponent. Such a slow convergence — if confirmed by modeling studies with higher resolution and additional physical processes — would provide some hope for longer-term, decadal-scale forecasts of the coupled ocean-atmosphere system.

## 4 Concluding remarks

We have presented a low-order model that incorporates several key physical ingredients of the coupled ocean-atmosphere system at mid-latitudes: baroclinic instability in a dry atmosphere, upper-ocean mass and heat transport in a shallow layer, an energy balance scheme incorporating radiative and heat fluxes, and a mechanical coupling mechanism between the atmosphere and the ocean. The latter coupling allows for the development of vigorous, albeit smooth oceanic gyres, cf. Figs. 4–5. The thermodynamic portion of the model represents a significant addition to previous model versions (73); (74). and is responsible for some of the most striking results reported herein.

This model contains variables that represent the dominant modes of variability of the coupled system, with variables for the atmosphere and for the ocean. Given the relatively large number of variables and processes involved, our model contains several

parameters that allow us to calibrate the intensity of the heat and momentum exchanges between the atmosphere and the ocean. The three most important ones are the wind stress coupling between the two fluid media, the meridional variation of gradient of the radiative flux coming from the sun, and the strength of the heat fluxes between the two.

Our bifurcation analysis and numerical integrations have revealed a rich variety of solutions and a diversity of regimes of model behavior; see Fig. 1. In particular, we found a coupled ocean–atmosphere mode for large values of the friction coupling and of the radiative-flux gradient . This coupled mode has a decadal time scale and it appears as a stable limit cycle for a range of parameter values, e.g. Wm at s; see green solution branch in Fig. 1a, the stable limit cycles in Fig. 3, as well as the purely periodic solutions plotted in Figs. 6 and 7.

Interestingly, the slow, decadal-scale solutions involve not just the oceanic variables, but some of the atmospheric ones as well, and correspond thus to true coupled modes. In fact, we found a subspace of dimension in the model’s 36-dimensional phase space that is invariant under the model flow and involves only slow motions. Even when the limit cycles in Fig. 3 lose their stability, they still seem to organize the phase space flow nearby into solutions that are dominated by slow motions; see the solutions that appear as small perturbations of a limit cycle in Figs. 6 and 7, as well as the upper panel of Fig. 8. All the solutions of this type are still dominated by a decadal and a bidecadal periodicity, as visually apparent in Fig. 8a and further documented by the SSA spectra in Fig. 9.

These results suggest that the slow limit cycles form the backbone of a strange attractor that is, at first, embedded in the slow, 17-dimensional subspace and only extends beyond it as the atmosphere gradually decouples from the ocean. In fact, already the results summarized in Fig. 7 show that, as the heat flux intensity between the two fluid media decreases, the solutions become more agitated, with the fast part of the variance increasing in strength.

We conjecture that this backbone corresponds to the least-unstable periodic orbit embedded in the chaotic attractor, as suggested by the discussion of Fig. 1 in (25). The term is justified here since this orbit can be clearly identified in our model’s singular spectrum as the main structure organizing its dynamics and a dominant source of its LFV.

We used the Kolmogorov-Sinai entropy to visualize in Fig. 12 that the perturbations of the slow limit cycles still contain large portions of relatively quiescent phase space flow, while other portions become more agitated. Next, we used the Lyapunov dimension to track the evolution of the strange attractor we suspect is embedded in the slow subspace, as parameters change. This was done for the three orbits plotted in Fig. 7 and for several additional ones, computed at the same and values.

For the slow limit cycle at WmK (blue curve in the figure), one obviously gets As decreases, increases, until it finally exceeds at In fact, when , i.e., when there is no heat flux at all between the two media, . Visually, the green and red objects in Fig. 7 look successively stranger, in full agreement with our Lyapunov dimension calculations.

We thus have substantial numerical evidence for the presence of a strange attractor in our coupled model, which is completely dominated by slow motions for sufficiently strong coupling between the atmosphere and ocean. As this coupling weakens, the dimension of the attractor increases and it supports more and more fast variance.

The slow variability associated with this attractor reflects the strong effects of the horizontal heat transport within the ocean on the atmospheric fields above. This horizontal transport intensifies, to first order, when the meridional temperature gradient in the ocean is large, and it slows down when the gradient is small. Concomitantly, a large meridional gradient in the ocean imposes large meridional gradients of heat flux and longwave radiation in the atmosphere; the latter, in turn, induce an increase of baroclinic instability and hence fast variability within the atmosphere, which enhance the heat transported toward the pole by the latter.

In the coupled model studied herein, these feedbacks between atmospheric and oceanic transports are activated by the mechanical coupling between the ocean and the atmosphere, via wind stress forcing — which induces the development of gyres — as well as by the vertical heat fluxes, which are proportional to . Both of these couplings, mechanical and thermal, are simplified here by modeling the temperature as a passive scalar. Our results on the effects of ocean–atmosphere coupling are plausible, given the fact that similar mechanisms have been found to act in GCMs (Delworth and Mann, 2000; L’Hévéder et al., 2014); still, it is well worth incorporating a more active coupling between ocean dynamics and temperature in future versions of the present model.

The periodicities found in our model are decadal and bidecadal, cf. Figs. 8a and Fig. 9. The former periodicity agrees roughly with the variability peak in the coupled modes of the slightly more highly resolved atmosphere–ocean model of Kravtsov et al. (2007, 2008). In the latter, the near-decadal peak was associated with alternations between two atmospheric regimes characterized by the subtropical jet’s latitude; this bimodality of the jet position and intensity does not seem to play a key role in the present model.

On the other hand, the NAO exhibits a peak in its variability at 7-8 years (Feliks et al., 2011, and references therein), rather than at 10–12 years, like here and in the coupled modes of Kravtsov et al. (2007, 2008). The mechanism that drives the 7-8-year peak in Feliks et al. (2011) is the purely oceanic gyre mode of (13); (28); (62) — which arises even for time-independent wind stress forcing — through this mode’s impact on the atmospheric jet that forms above the SST front associated with the gyre mode, cf. Feliks et al. (2004); Feliks et al. (2007); (4).

To clarify further the nature of the decadal and multidecadal variability in the present coupled model — and its relation, if any, with the LFV in and around the North Atlantic — will require two things. First, the spatio-temporal aspects of the slow coupled mode will have to be examined using multi-channel SSA (25), rather than just the single-channel version applied herein. The results of such an analysis can then be compared more closely with the observed spatio-temporal patterns. Second, additional modes, zonal as well as meridional, will have to be added within both the oceanic and the atmospheric component of the low-order model in order to clarify the robustness of the LFV to a more detailed spatial representation.

The stability of the coupled model’s phase space flow was measured by its leading Lyapunov exponents (Figs. 10, 13 and 14), as well as by its local Kolmogorov-Sinai entropy (Fig. 12). We found that these stability measures are highly dependent on the parameter values, in particular on the coupling parameter . For small values of , the transport of heat within the ocean is relatively small and the coupled system’s instability is essentially driven by the atmosphere; hence the Lyapunov exponents are quite large. Once attains a certain threshold — which depends in turn on the values of and of — the values of the Lyapunov exponents are drastically reduced, cf. Fig. 13, while the oceanic transport intensifies.

The drop in the coupled model’s short-term error growth is obviously associated with an increase in predictability; this increase, in turn, emphasizes the ocean’s potential importance in modulating our ability to predict both components of the coupled system at short lead times. Note that this oceanic effect was not present in the previous, 24-equation version of the model (74).

The main difference between the two model versions is the absence of a temperature equation in the ocean and of an energy balance between the ocean and the atmosphere in the earlier version: in (74), the thermal energy source is provided through a restoring force toward climatological radiative equilibrium within the atmosphere. In the earlier version there was, therefore, no heat transport within the ocean, and in turn no thermal impact on the atmosphere. As the present analysis shows, this tranfer of energy plays a crucial role in the predictability of the coupled system.

At long lead times, the coupled model’s predictability is essentially dominated by the role of the LFV (see again Fig. 8) and of the hypothetical slow attractor, cf. Fig. 12. When the slow attractor is present long-lead forecasts can be fairly accurate even for the atmospheric variables; see Fig. 15. It is, therefore, worth investigating whether such slow, coupled modes are indeed found within the real coupled ocean-atmosphere system or, at least, within a much more detailed coupled model.

Clearly, our relatively low-order, 36-dimensional coupled model has produced a substantial number of interesting and stimulating results. Several of its limitations have already been mentioned and we summarize the main ones here.

First, we have been working with a closed, rectangular ocean basin, somewhat similar to those in (28); (62). The presence of oscillations between the cyclonic and anticyclonic gyres in the upper ocean for such a configuration is by now well known, although its robust presence in much more realistic intermediate models has also been demonstrated ((13), and references therein).

Second, the emphasis on the wind-driven, double-gyre circulation neglects the importance of the buoyancy-driven, overturning circulation (13); (58); (63) for decadal and multi-decadal climate variability (6); (20). Some aspects of both types of oceanic circulation were included, for instance, in Chen and Ghil (1996) and in Kravtsov et al. (2007, 2008). A related caveat is the passive scalar character of the temperature, which does not feed back on the ocean dynamics.

Among the missing processes in our model formulation is the hydrological cycle, which can affect substantially the atmospheric variability (42), as well as the heat fluxes between the atmosphere and the ocean (4); Chen and Ghil (1996). To address the robustness of some of the key results obtained herein — including the presence of slow coupled modes and of their impact on predictability — will require removing some of these limitations and expanding further the methodology used in the present work.

## 5 Acknowledgments

This work originated in discussions between MG and SV during the Mathematics for the Fluid Earth Programme held at the Isaac Newton Institute in Cambridge, UK, in Fall 2013. We have also benefited from discussions with C. Nicolis and H. Goosse. This work is partially supported by the Belgian Federal Science Policy Office under contracts SD/CA/04A and BR/12/A2/STOCHCLIM (SV, JD, LDC), by the U.S. National Science Foundation under grant OCE-1243175 (MG), and by the U.S. Office of Naval Research under MURI grant N00014-12-1-0911 (MG).

## Appendix A Nondimensional equations: Dynamics

The equations (1)–(3) are nondimensionalized by scaling horizontal distances by , time by , the vertical velocity by and the atmospheric and oceanic streamfunctions and by . The parameters are also rescaled as , , , and .

The atmospheric and oceanic fields are expanded in Fourier series over the domain , where is the aspect ratio between the meridional and the zonal extents of the domain, . The values of the parameters just discussed above that are held fixed in the rest of the paper are km, s, , , , , and ; the value of the Coriolis parameter corresponds to the domain’s axis of N–S symmetry at N.

## Appendix B Linearization around a climatological temperature

We assume here that

(10a) | |||

(10b) |

where and are climatological reference temperatures, and . Then

(11a) | |||

(11b) |

We also let

(12a) | |||

(12b) |

with and the averaged shortwave radiative forcings.

Neglecting the higher-order terms in in Eq. (11b) and separating the reference and the perturbation — i.e., the zeroth-order terms in the expansion from the first-order ones — leads to the following temperature equations for the ocean:

(13a) | |||

(13b) |

Applying the same procedure to the atmospheric temperatures yields:

(14a) | |||

(14b) |

It is interesting to note that the equations (13a, 14a) for the climatological references are independent of the perturbations. This implies that stationary solutions can be readily found by solving

(15a) | |||

(15b) |

Summing the two equations and substituting the expression obtained for into the first one leads to

(16) |

a quartic equation for the mean ocean temperature . Solving this algebraic equation analytically by using Mathematica ((77)) leads to a unique real and positive stationary solution for , a typical ratio for the shortwave radiation absorbed by the atmosphere vs. that absorbed by the ocean. Moreover, this stationary solution is stable for a wide range of values of .

The stationary solutions so obtained, however, are not really close to the observed values for the coupled ocean–atmosphere system. In order to have our model operate in a more realistic domain of phase-parameter space, we set K and K, as in (3). This choice will only affect the values of the parameters in front of the linearized radiative terms of Eqs. (13) and (14).

## Appendix C Derivation of the truncated equations

In this appendix, we describe how the low-order model is derived, by projecting the equations presented in section 2 onto a suitable orthogonal basis.

For the atmosphere, we keep the same set of modes as in (53); (74),

(17) |

while the imposed radiative fluxes are given by and One furthermore assumes that , where is the fraction of shortwave radiation absorbed by the atmosphere. Note that two values of this fraction will be used, and 1/3, both of which lie in a range that is typical for Earth’s atmosphere.

For the ocean, we retain the following set of modes:

(18) |

This set contains two additional modes in the latitudinal direction, as compared to Veronis (1963); see also ((28); (59)). This addition will allow us to get a more realistic temperature profile within the ocean.

We assume hereafter that one can project and truncate the atmospheric and oceanic streamfunction fields and , as well as the corresponding temperature fields and , and the vertical velocity according to

(19a) | |||

(19b) | |||

(19c) |

where and , with and the streamfunctions in the upper and lower layer of the atmosphere.

The projection of the dynamical equations (1)–(3) on the modes kept in Eqs. (C) and (C) has already been discussed in detail in (53); (74), and it will not be repeated here. We describe, however, in detail the procedure for projecting the temperature equations.

Non-dimensionalizing the baroclinic and barotropic streamfunction fields, and , and the parameters of the truncated, low-order model, such that , , , and , one gets

(22) | |||||