Multiple Equilibria and Oscillatory Modes in a Mid-latitude Ocean-forced Atmospheric Model

Atmospheric response to a mid-latitude sea surface temperature (SST) front is studied, while emphasizing low-frequency modes induced by the presence of such a front. An idealized atmospheric quasi-geostrophic (QG) model is forced by the SST field of an idealized oceanic QG model. First, the equilibria of the oceanic model and the associated SST fronts are computed. Next, these equilibria are used to force the atmospheric model and compute its equilib-ria when varying the strength of the oceanic forcing. Low-frequency modes of atmospheric variability are identified and associated with successive Hopf bifurcations. The origin of these Hopf bifurcations is studied in detail, and connected to barotropic instability. Finally, a link is established between the model's time integrations and the previously obtained equilibria.


Introduction
Over the past decades, many studies have addressed the statistical and synoptic identification and description, as well as the physical explanation of atmospheric weather regimes, as characterized by recurrent and persistent circulation patterns.Using statistical methods, it has been demonstrated that the two opposite phases of the North Atlantic Oscillation -blocked and zonal, respectively -are the stablest and thus most commonly accepted weather regimes in the Atlantic sector (Cheng and Wallace, 1993;Kimoto and Ghil, 1993;Corti et al., 1999;Smyth et al., 1999;Cassou et al., 2004).If, however, a prevailing consensus on the observed weather regimes has been established, there is still a debate on the explanation of their occurrence.
When looking at a map of atmospheric variability, we observe that its Northern Hemisphere maxima are located over the Atlantic and Pacific Oceans, thus suggesting a possible link between the sea surface temperature (SST) field and the atmospheric variability.The response of the atmosphere to mid-latitude SST anomalies has been a contentious issue over the last few decades (Kushnir et al., 2002).Recent work (Feliks et al., 2004(Feliks et al., , 2007;;Minobe et al., 2008;Small et al., 2008;Hogg et al., 2009;Chelton and Xie, 2010) has shown that oceanic SST fronts play a key role in this response.Sweet et al. (1981) and Businger and Shaw (1984) already noted that an SST front is responsible for the unequal heating of the atmosphere's lower layers on either side of the front.This differential heating is transmitted to the troposphere via several processes (Lee and Kim, 2003;Nakamura et al., 2008;Deremble et al., 2012).As a result, the growth of baroclinic eddies and the position of the jet stream are strongly affected by the position and strength of an SST front.
In order to explain the dynamic origin and maintenance of multiple weather regimes, Charney and Devore (1979) suggested to compute the steady states of an idealized atmospheric model.They found a clear link between blocked and zonal weather regimes and the multiple equilibria of their model.Similar results have been obtained with more complex and realistic models (Reinhold and Pierrehumbert, Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union. 1982; Legras and Ghil, 1985;Crommelin, 2003;Sempf et al., 2007;Deremble et al., 2009).These studies mostly focused on the steady states and low-frequency oscillatory modes of spectrally truncated atmospheric models.
The mid-latitude oceanic circulation has also been studied by applying the methods of numerical bifurcation theory, leading to the identification of multiple steady states, as well as of low-frequency oscillatory modes.Indeed, it has been shown that many aspects of the interannual variability of the double-gyre circulation in the mid-latitude ocean basins can be explained in terms of steady states and of irregular, chaotic oscillations of the oceans' evolution equations, forced by time-constant, climatological winds (see Dijkstra, 2005 or Dijkstra and Ghil, 2005 for a review).The gyre modes that arise in this kind of study provide a possible explanation for oceanic low-frequency modes (Jiang et al., 1995;Speich et al., 1995;Simonnet and Dijkstra, 2002).
The purpose of the present paper is to apply the fruitful approach of successive bifurcations to the study of the atmosphere's intrinsic variability in the presence of oceanic SST fronts.Both our model's atmospheric and oceanic components are modeled by the quasi-geostrophic (QG) barotropic equation, in a periodic β-channel and an idealized rectangular basin, respectively.The SST is used here as the forcing that drives the free troposphere via the vertical velocities that are calculated analytically at the top of the marine atmospheric boundary layer (MABL).We examine first how the atmosphere reacts to the SST associated with the distinct oceanic equilibria.Next, we isolate the atmospheric modes of low-frequency variability (LFV) to which different SST fields give rise.To do so, we search for the Hopf bifurcations resulting in these modes.
Our main focus is the atmosphere's intrinsic variability, and we leave the fully coupled model for future study.The study of LFV in more or less idealized, fully coupled oceanatmosphere models (Kravtsov and Robertson, 2002;Berloff et al., 2007a;Kravtsov et al., 2007a,b) has already started and will not be discussed here any further.
We emphasize in our model of ocean-atmosphere interaction the pressure gradient anomaly created by variations of the air temperature in the MABL, known as the thermodynamic pressure gradient, and do not take into account explicitly the effect of the vertical mixing on the height of the boundary layer (e.g., Wallace et al., 1989;O'Neill et al., 2005;Samelson et al., 2006).Lindzen and Nigam (1987) have parametrized this process in a tropical boundary layer, and several studies (e.g., Battisti et al., 1999;Bellon et al., 2003;Chen et al., 2003) have applied their parametrization in the tropics.Feliks et al. (2004) (FGS hereafter) bypassed the Lindzen and Nigam (1987) calculation in the extratropics and showed that the thermal component of the vertical velocity at the top of the MABL is associated with the Laplacian of the SST field.The FGS parametrization has been applied over the North Atlantic basin (Feliks et al., 2011), where it reproduces correctly major aspects of the atmospheric circulation.Brachet et al. (2012) have further verified this parametrization in a general circulation model (GCM) with zooming over the Gulf Stream region and found it to represent the dominant effects of ocean-atmosphere interaction there.
Of course, air-sea interactions involve much more complex phenomena than the FGS parametrization chosen here (see Small et al., 2008 for an extensive review).The model behavior is thus not expected to reproduce the observations in great detail.Our intent is merely to extract from this idealized model the minimal ingredients that permit one to highlight some generic features of ocean-atmosphere interaction.
The paper is organized as follows.In Sect.2, we give a brief description of the model.In Sect.3, we compute the ocean's multiple equilibria as a function of the intensity of the time-constant wind field.In Sect.4, we study the atmosphere's multiple equilibria as a function of the intensity of the SST front computed in the oceanic model; this study is extended to time-dependent solutions of the atmospheric model in Sect. 5.The study of the fully coupled model is left for a subsequent paper.Results are summarized and discussed in Sect.6.

Model description
We construct an ocean-atmosphere model with an MABL that helps parametrize the momentum fluxes between the ocean and atmosphere.The model's atmospheric and oceanic components are each governed by a barotropic QG equation.Our model's free atmosphere and MABL follow FGS, while the oceanic component follows Simonnet and Dijkstra (2002) and Simonnet et al. (2005).We recall that the effects of the atmosphere on the ocean, as well as those of the ocean on the atmosphere, are studied independently.Full coupling between the two is deferred for a further study.
A schematic model diagram is shown in Fig. 1.There are three layers: the oceanic layer of height H oc , and the two atmospheric layers, of total height H at , with an Ekman layer of height H E and the free troposphere of height H at −H E .In this figure, we plot the fluxes between the three layers and some of the source terms of these fluxes.The term that forces the atmospheric vorticity is composed of two parts: a classical dissipation term, γ ∇ 2 ψ at , and a forcing term proportional to the Laplacian of the SST field, α∇ 2 T oc .There is no external forcing of the atmospheric layer, such as a differential poleto-equator heating.
In contrast, the ocean is forced by a time-constant wind stress term σ F, which mimics the climatological surface winds at midlatitudes, and subject to bottom friction with a constant relaxation time τ ek .In addition, the SST is relaxed to an equilibrium profile chosen to resemble the temperature profile of extratropical oceanic basins; this SST forcing term is not shown in the figure.The z-axis points upward, while the x-and y-axes point east-and northward, respectively.

The free-atmosphere model
The evolution of the atmospheric potential vorticity (PV), q at = ∇ 2 ψ at , is formulated by using the QG approximation on the β-plane, β being the meridional gradient of the Coriolis parameter, and ψ at the atmospheric stream function (see Fig. 1).Note that the stream function is linked to the pressure since P (H E ) = ρ 0 f 0 ψ at .This PV equation is integrated for a single layer whose height is H at − H E and with a rigid lid, i.e., with zero vertical velocity at the top of the layer and an infinite deformation radius.In non-dimensional variables, we obtain with J (•, •) the Jacobian operator and R at the atmospheric Reynolds number; where L r and U r are a characteristic length and velocity, respectively, and ν at is the horizontal eddy viscosity.
The last two terms of Eq. ( 1) correspond to the friction and the forcing, respectively (see also Appendix B).The seaair forcing parameter α and the relaxation constant γ correspond to the non-dimensional proportionality coefficients in Eq. (B8): To nondimensionalize the equation for γ , we used the Froude number F r = U r /(gH E ) 1/2 .

The barotropic ocean model
Similarly, we write the evolution equation for the oceanic PV in non-dimensional form: with ψ oc the oceanic stream function.The oceanic PV q oc is computed while taking a free upper surface, instead of the rigid lid taken for the atmospheric part: where R d is the deformation radius of the barotropic ocean, while the Reynolds number R oc is defined using the eddy viscosity ν oc of the ocean.The third term on the right-hand side of Eq. ( 4) is the bottom friction, with a time constant that is equal to τ ek .The forcing F corresponds to the non-dimensional wind stress curl and it is written as with y = 0 at the south wall of the basin and y = 1 at the north wall, in non-dimensional units.This forcing corresponds to polar and tropical easterlies near the ocean basin's northern and southern walls, and to westerlies near its midlatitude symmetry axis.

The sea surface temperature (SST) field
The oceanic temperature T oc is relaxed towards an imposed profile T , and also advected by the oceanic velocity field as a passive tracer.Its full evolution equation is given by with k H the thermal diffusivity, and χ an inverse time constant for the thermal relaxation.The profile T = T (y) toward which we relax the temperature is given by with T the north-to-south temperature difference in the basin and η a constant that calibrates the steepness of the relaxation profile.

Parameter choices
We list all the parameters used for the basic, or standard, model solution in Table 1.Some parameters -like f 0 , ρ at , ρ oc -are kept fixed, while other parameters -such as the eddy viscosities ν at and ν oc or the bottom friction parameter (1/τ ek ) -affect the behavior of the fluid: depending on the value chosen, the fluid will be more or less turbulent.The standard values chosen here (cf.Table 1) are dictated by the horizontal resolution of 30 km.The thermal relaxation constant χ is set so as to keep the relaxation term small with respect to the advection term in the basin's western part.We will focus mainly on the wind intensity σ and the height of the layers H E , H at , and H oc and use the non-dimensional version of the parameters.As noted in the Introduction, these heights can vary quite widely depending on the atmospheric or sea surface conditions.We explore here the parameter plane (σ, α), where σ is the strength of the climatological wind-stress forcing.By varying this parameter, we verify in Sect. 3 the well-known results on multiple equilibria of the wind-driven ocean circulation, and thus calibrate our model.In Sects. 4 and 5, we vary the forcing α of the SST field on the atmosphere, in order to see how the latter reacts to a change in the strength and zonal extent of the SST front.

Numerical scheme and discretization
For the set of partial differential equations that govern our model's atmosphere, ocean, and SST field -namely Eqs.(1), ( 5) and (7) -the boundary conditions are -periodic boundary conditions at the eastern and western ends of the β-channel, and free-slip conditions along its northern and southern walls; -free-slip conditions for the oceanic stream function on all basin walls; and -no-flux boundary conditions for the temperature, i.e., its derivative is zero perpendicular to each boundary.
We discretize the variables on a A-grid (Arakawa, 1966) with the horizontal mesh size fixed at 30 km for the oceans as well as the atmosphere.This is an acceptable resolution to observe phenomena associated with frontal dynamics in the atmosphere (FGS; Minobe et al., 2008).
We use an LU factorization to invert the discrete Laplacian operator that converts the stream function ψ into the PV q.For this purpose, we use the UMFPACK library that provides an iterative solver of LU type (Davis, 2004).For the time integration, we use a third-order Adams-Bashforth scheme, which has good stability properties when the integration step is sufficiently small (Durran, 1999); we fix this step here to 15 minutes.

Multiple equilibria of the oceanic circulation
Many authors have studied multiple equilibria of the winddriven circulation using various models and diverse configurations (e.g., Veronis (1963), Cessi and Ierley (1995), Jiang et al. (1995), Speich et al. (1995), Primeau (1998) and Simonnet et al. (2005)).We present here the corresponding results for our model as a basis for the rest of our study.The barotropic ocean in this section is driven by the timeconstant, idealized wind profile of Eq. ( 6), and we study the solutions as the intensity σ of the wind stress increases.
Figure 2 shows the bifurcation diagram obtained by varying σ , all other parameters in Table 1 being kept constant.This figure has been obtained using the pseudo-arclength continuation method of Keller (1977), as introduced into atmospheric studies by Legras and Ghil (1985) and into oceanographic ones by Speich et al. (1995) (see Dijkstra (2005) for numerical details).We plot on the ordinate the energy difference E oc between the subtropical gyre and the subpolar gyre: E = ψ>0 ψ 2 − ψ<0 ψ 2 .This diagram is particularly robust, as it resembles the ones obtained in all the above-mentioned studies of the double-gyre, wind-driven circulation.
Among the robust elements -appearing in all the studies that use perfectly symmetric, QG models -we note the presence of an antisymmetric solution for which the subtropical and the subpolar gyre have the same strength, for all σ -values.In shallow-water models, the mirror symmetry with respect to the east-west axis of the basin is broken and the pitchfork bifurcation takes the perturbed form shown by Jiang et al. (1995).Still, for a wide range of σ -values, one finds a nearly antisymmetric flow pattern.respect to the East-West axis of the basin is broken and the pitchfork bifurcation takes the perturbed form shown by Jiang et al. (1995).Still, for a wide range of σ-values, one finds a nearly antisym-185 metric flow pattern.
For the value of the oceanic Reynolds number Roc = 5 × 10 3 in Table 1, there are two successive pitchfork bifurcations, which give respectively two branches of asymmetric equilibria, symmetric with respect to each other.The bifurcation at P1 is supercritical, as the antisymmetric branch loses its stability to the two asymmetric ones, while the one at P2 is subcritical, as the antisymmetric 190 branch regains its stability, while the two asymmetric ones that grow out of it are unstable.
In Fig. 3, we plot an asymmetric solution for a value of σ just below the first Hopf bifurcation, H ′ 1 (∆Eoc > 0).The asymmetry of the solution in panel (a) is apparent from the fact that the subtropical 9 Fig. 2. Bifurcation diagram of the oceanic circulation.The energy difference E oc between the subpolar and the subtropical gyre is plotted as a function of the wind intensity σ in nondimensional units; all other model parameters equal their values in Table 1.The square P 1 marks a supercritical pitchfork bifurcation and the square P 2 a subcritical pitchfork bifurcation, the equilibrium branch between the two being unstable.Given the large number of unstable branches in the atmospheric model (see Fig. 5 in Sect. 4 below), we prefer for clarity's sake not to use the customary dashes for the unstable branches.The circles H 1 , H 2 , H 3 , H 4 and H 1 , H 2 , H 3 , H 4 mark the position of successive supercritical Hopf bifurcation points, at which either equilibrium branch loses the stability of yet another complex conjugate eigenvalue (see text for details).
For the value of the oceanic Reynolds number R oc = 5 × 10 3 in Table 1, there are two successive pitchfork bifurcations, which give respectively two branches of asymmetric equilibria, symmetric with respect to each other.The bifurcation at P 1 is supercritical, as the antisymmetric branch loses its stability to the two asymmetric ones, while the one at P 2 is subcritical, as the antisymmetric branch regains its stability, while the two asymmetric ones that grow out of it are unstable.
In Fig. 3, we plot an asymmetric solution for a value of σ just below the first Hopf bifurcation, H 1 ( E oc > 0).The asymmetry of the solution in Fig. 3a is apparent from the fact that the subtropical gyre is stronger than the subpolar one.At the same time, the spatial extent of the subpolar gyre is greater than that of the subtropical gyre; see also Chang et al. (2001).
Associated to the oceanic circulation in Fig. 3a, the SST field in Fig. 3b shows that the front is distorted as it follows the contours of the stream function and is thus deflected toward the south.In the eastern part of the basin, the temperature profile closely resembles the imposed forcing, since the ocean dynamics is not very active in the so-called Sverdrup interior, away from the intense recirculation near the western boundary.

B. Deremble et al.: Multiple equilibria and oscillatory modes in a mid-latitude atmosphere model
In the eastern part of the basin, the temperature profile closely resembles the imposed forcing, since the ocean dynamics is not very active in the so-called Sverdrup interior, away from the intense recirculation near the western boundary.On the asymmetric branches emerging from the first pitchfork bifurcation, we note the successive appearance of four Hopf bifurcations; they correspond to oscillatory instabilities with periodicities that equal, respectively, 800, 200, 300 and 1200 days.The spatio-temporal patterns of these instabilities can be illustrated by plotting the real and imaginary parts of the eigenvector associated with each of them, as shown in Fig. 4 for the first of them, at H 1 .
The two plots, in Fig. 4, correspond to two distinct phases of the oscillatory mode.This mode is concentrated in the basin's western part and corresponds to an alternative strengthening and weakening of each gyre; it has been labeled the gyre mode, and has been extensively studied (see Simonnet and Dijkstra, 2002).This oscillatory mode also corresponds to a displacement of the front around its equilibrium position; its period here is 800 days or roughly 2.2 yr.
Such modes are of great interest because they are likely to be responsible for the low-frequency, sub-and interannual variability of the mid-latitude ocean.Chang et al. (2001), Nadiga and Luce (2001) and Simonnet et al. (2005) showed that a homoclinic bifurcation appears at a slightly higher value of σ .This global bifurcation in turn generates very low-frequency, interdecadal oscillations.At this point, the system is able to switch from one asymmetric circulation to the opposite one.Simonnet and Dijkstra (2002) were interested in describing more precisely the Hopf bifurcation associated with the gyre mode.These authors found that the eigenvalue crossing the imaginary axis, and responsible therefore for the instability, arises from the merging between the eigenvalue associated with the first pitchfork bifurcation, at P 1 , and another eigenvalue associated with a saddle-node bifurcation on the main, antisymmetric branch at a larger value of σ .The former has an asymmetric eigenfunction that causes the symmetry breaking; the latter is a shear instability.At the merging point, the imaginary part of the eigenvalue is small and the appearance of four Hopf bifurcations; they correspond to oscillatory instabilities with periodicities that equal, respectively, 800, 200, 300 and 1200 days.The spatio-temporal patterns of these instabilities can be illustrated by plotting the real and imaginary parts of the eigenvector associated with each of them, as shown in Fig. 4 for the first of them, at H ′ 1 .

205
The two plots, in panels (a) and (b), correspond to two distinct phases of the oscillatory mode.
This mode is concentrated in the basin's western part and corresponds to an alternative strengthening and weakening of each gyre; it has been labeled the gyre mode, and has been extensively studied (see Simonnet and Dijkstra, 2002).This oscillatory mode also corresponds to a displacement of the front around its equilibrium position; its period here is 800 days or roughly 2.2 years.

210
Such modes are of great interest because they are likely to be responsible for the low-frequency, sub-and interannual variability of the mid-latitude ocean.Chang et al. (2001), Nadiga and Luce (2001) and Simonnet et al. (2005) showed that a homoclinic bifurcation appears at a slightly higher value of σ.This global bifurcation in turn generates very low-frequency, interdecadal oscillations.
At this point, the system is able to switch from one asymmetric circulation to the opposite one.
215 Simonnet and Dijkstra (2002) were interested in describing more precisely the Hopf bifurcation associated with the gyre mode.These authors found that the eigenvalue crossing the imaginary axis, and responsible therefore for the instability, arises from the merging between the eigenvalue associated with the first pitchfork bifurcation, at P 1 , and another eigenvalue associated with a saddlenode bifurcation on the main, antisymmetric branch at a larger value of σ.The former has an 220 asymmetric eigenfunction that causes the symmetry breaking, the latter is a shear instability.At the merging point, the imaginary part of the eigenvalue is small and the associated mode is of lowfrequency type.Hence the gyre modes (i) have low frequency; and (ii) are quite distinct from the 11 associated mode is of low-frequency type.Hence, the gyre modes (i) have low frequency; and (ii) are quite distinct from the classical basin Rossby modes of Sheremet et al. (1997).
As the domain or the basic-state flow is deformed, cascades of symmetry-breaking transitions may occur.The beginning of such a cascade is indeed observed here, and it is described in Appendix A. For example, in the context of jetstream penetration into an oceanic basin, this cascade corresponds to the occurrence of stationary Rossby waves with increasing wavenumber and nontrivial dispersion relations (see also Primeau (1998)).More generally, Simonnet (2008) has shown that the bifurcations associated with large-scale patterns, as observed here, arise from the discrete spectral behavior of the 2-D Euler equations.In particular, these bifurcations are not dependent on a particular eddy-viscosity parametrization, as it is often claimed.Simonnet (2005) argued that the link between the symmetry breaking -and therefore the presence of multiple equilibria -and the emergence of low-frequency modes is part of a natural and logical response to increasing PV flux from the wind stress: the ocean circulation is first distorted from its original symmetry so as to draw more PV from the forcing in steady state.Next, the only way to react to a further increase in PV input is by becoming unstable (see also Ghil et al. (2008) for a broad outlook on this interpretation of LFV generation in geophysical flows).Our study of atmospheric LFV in the next section shows that similar phenomena occur there, too, and will help clarify the reasons for the instability's low frequency.

Multiple atmospheric equilibria
In this section, we vary the parameter α that governs the forcing of the SST field on the free atmosphere.We start our investigation with the multiple oceanic equilibria that have been found in the previous section and study the associated atmospheric equilibria as the parameter α varies.
According to Eq. ( 3), we see that α is proportional to the square of the height H 2 E of the Ekman layer, α ∝ H 2 E .For reference purposes, a height H E = 300 m corresponds to α = 0.04 (see Table 1), while H E = 1000 m yields a value of α = 0.5.The height of the boundary layer is also difficult to estimate accurately over large areas, according to the studies reviewed in the introduction.

Short SST front
For each antisymmetric equilibrium of the oceanic circulation, we now vary the parameter α.We begin with an oceanic front of modest longitudinal extent -about 300 km -which is obtained for σ 0.6.This value corresponds to an SST front similar to that plotted in Fig. 3c.Referring to the bifurcation diagram in Fig. 2, we see that such an antisymmetric front is no longer stable at this parameter value.In fact, at σ = 0.55, we are already past the first pitchfork bifurcation, P 1 , that destabilizes the antisymmetric equilibrium.
For σ > 0.55 , there are two stable equilibria that correspond both to asymmetric ocean circulations, and are mirror images of each other.The strength of such an SST front is roughly 3 K/10 km.This value is comparable to the observed Gulf Stream or Kuroshio front -cf.FGS, Feliks et al. (2007), Feliks et al. (2011) and Brachet et al. (2012) -although the impact of eddies does affect this value in observed western boundary currents.
The bifurcation diagram for the atmospheric equilibria with respect to α is shown in Fig. 5, for σ = 0.65.We plot, as in Fig. 2 for the oceanic case, the energy difference E at between the cyclonic and the anticyclonic circulation.As for the ocean's double-gyre circulation, there exists an antisymmetric equilibrium for all α-values.This antisymmetric equilibrium is almost always unstable since, even at very  The bifurcation diagram for the atmospheric equilibria with respect to α is shown in Fig. 5, for σ = 0.65.We plot, as in Fig. 2 for the oceanic case, the energy difference ∆Eat between the cyclonic and the anticyclonic circulation.As for the ocean's double-gyre circulation, there exists an antisymmetric equilibrium for all α-values.This antisymmetric equilibrium is almost always unstable since, even at very small α-values, there appear several Rossby-wave instabilities, with a periodicity of about 7 days, along with the associated harmonics.These high-frequency modes are not plotted in Fig. 5, because we are only interested here in low-frequency phenomena.Thus, in contrast to the oceanic equilibria of Fig. 2, the atmospheric equilibria in Fig. 5 are almost always unstable.
Two pitchforks bifurcations, at α = 0.07 and α = 0.14, arise next.It follows that multiple equilibria of the atmospheric circulation are not only generated by its interaction with actual, "mechanical" 13 Fig. 5. Bifurcation diagram for the atmospheric equilibria at σ = 0.65 .The asymmetry measure E at is computed in the same manner as in Fig. 2 for the ocean but as a function of the parameter α.
All branches, whether stable or not, are shown as solid lines, for legibility.On the steady-state branches, we only plot -with circles -the Hopf bifurcations for which the period exceeds 50 days.P 1 is a supercritical pitchfork bifurcation, and P 2 is a subcritical one.
The shape of the SST forcing used here is similar to Fig. 3c.
small α-values, there appear several Rossby-wave instabilities, with a periodicity of about 7 days, along with the associated harmonics.These high-frequency modes are not plotted in Fig. 5, because we are only interested here in lowfrequency phenomena.Thus, in contrast to the oceanic equilibria of Fig. 2, the atmospheric equilibria in Fig. 5 are almost always unstable.Two pitchforks bifurcations, at α = 0.07 and α = 0.14, arise next.It follows that multiple equilibria of the atmospheric circulation are not only generated by its interaction with actual, "mechanical" topography, as in Charney and Devore (1979) or Legras and Ghil (1985), but also with the "thermal" topography of SST fronts, through the vertical velocities and PV fluxes generated by the latter.These two bifurcations correspond to the same eigenvalue that becomes positive between the two values of α cited above.
Beyond the second pitchfork bifurcation P 2 , several asymmetric equilibria exist, up to high α-values.These branches, however, are not connected to the main branch via a pitchfork bifurcation; they are, therefore, quite difficult to detect and locate.Their presence lets us suspect that the bifurcation diagram in Fig. 5 might be incomplete, since the continuation method used herein does not allow one to find all the disconnected branches.
As α increases, the asymmetric component of the steadystate circulations grows, too, as is the case for the oceanic equilibria, when increasing σ .According to Fig. 5, there exist between one and four pairs of asymmetric equilibria for α > 0.07.
In Fig. 6, we display the spatial pattern of an antisymmetric and an asymmetric equilibrium for α = 0.12.The two solutions have several common characteristics.They are both dominated by a zonal or nearly zonal eastward jet, close to the channel's east-west axis of north-south antisymmetry.Westerly winds prevail along both the northern and southern walls of the β-channel.In the asymmetric solution, the jet is slightly wavy, and this undulation extends well past the SST front, away from the western portion of the channel.In both cases, the maximum jet strength is 8 m s −1 .This anomaly is of reasonable size and must be compared to the observed jet stream in the Atlantic, whose strength is 30-50 m s −1 (Peixoto and Oort, 1992).Indeed, as emphasized by FGS and Feliks et al. (2007), the jet that is due to the oceanic SST front is merely superimposed on the one caused by the classical eddy convergence of Lorenz (1967).

Characteristics of an unstable oscillatory mode
On the asymmetric branches originating from P 1 in Fig. 5, we note the presence of two Hopf bifurcations: H 1 and H 2 .As α increases, H 1 is a transition from an unstable fixed point -unstable, as previously mentioned, because of highfrequency Rossby waves arising at lower α -to a limit cycle; H 2 is the reverse-type transition, with a limit cycle shrinking to give a fixed point that is not stable either.In this sense, both Hopf bifurcations are supercritical: the first one in the conventional left-to-right direction of increasing bifurcation parameter, and the second one in the reverse direction.
These two bifurcations are associated with a pair of complex conjugate eigenvalues crossing the imaginary axis in one direction (from negative to positive) and then in the opposite direction (from positive to negative).The period of the limit cycles at the bifurcation point is of 145 days for H 1 and of 78 days for H 2 .We have not found any low-frequency Hopf bifurcations on the secondary, disconnected branches.
Figure 7 shows the spatio-temporal pattern associated with these two Hopf bifurcations; using the real and imaginary parts of the unstable eigenvector, one can build the different phases of the oscillation.As seen from the figure, the instability corresponds, in both cases, to a standing oscillation, in shape as well as intensity, of the atmospheric jet around its equilibrium position.In both cases, this oscillation is associated with a wavenumber-1 pattern, which is highly localized near the channel's symmetry axis.Therefore, we call this instability the jet mode.The spatial pattern of the two oscillatory instabilities is quite similar, which is not surprising, given the fact that the two are associated with the same pair of eigenvalues that crosses the imaginary axis, first to the right and then back to the left of it.

Origin of the instability
Following Simonnet and Dijkstra (2002) in their study of the oceanic circulation's gyre mode, we wonder next whether case for the oceanic equilibria, when increasing σ.According to Fig. 5 there exist between one and four pairs of asymmetric equilibria for α > 0.07.In Fig. 6, we display the spatial pattern of an antisymmetric and an asymmetric equilibrium for α = 0.12.The two solutions have several common characteristics.They are both dominated by a zonal or nearly zonal eastward jet, close to the channel's East-West axis of North-South antisym-285 metry.Westerly winds prevail along both the northern and southern walls of the β-channel.In the asymmetric solution, the jet is slightly wavy, and this undulation extends well past the SST front, away from the western portion of the channel.In both cases, the maximum jet strength is 8 m s −1 .This anomaly is of reasonable size and must be compared to the observed jet stream in the Atlantic, whose strength is 30-50 m s −1 (Peixoto and Oort, 1992).Indeed, as emphasized by FGS and Feliks We leave for further study the tracking of the complex conjugate pair that arises at M 1 , to verify there is not also a relationship between the occurrence of this atmospheric jet mode and the pitchfork bifurcation that precedes it.To clarify this issue, we plot in Fig. 8 the leading real eigenvalue λ 0 obtained when following one of the two asymmetric branches, from P 1 on.Tracking λ 0 is relatively difficult because of the density of complex conjugate eigenvalues whose real part is negative and that are close to the imaginary axis.Figure 8, however, demonstrates that there is no ambiguity in the monitoring of λ 0 : it becomes negative as soon as α passes its value at P 1 , continues to decrease with increasing α, and finally merges with another real eigenvalue (not shown) to form a complex conjugate pair.The merging takes place at the point M 1 in Fig. 8.
We plot in Fig. 9 the evolution of the spatial pattern of the eigenvector associated with this eigenvalue, from the pitchfork bifurcation P 1 until the merging M 1 .Figure 9a shows a purely symmetric eigenvector, associated in fact, at this point, with the pitchfork bifurcation.As α increases, the asymmetric component of the eigenvector grows, too.
We leave for further study the tracking of the complex conjugate pair that arises at M 1 , to verify that it is indeed responsible for the Hopf bifurcation at α = 0.09.The comparison of Fig. 9c with Fig. 7 is not conclusive.A clue that seems, however, to confirm this hypothesis is the smallness of the imaginary part of this eigenvector, at each of the two Hopf bifurcations, in comparison to the real part: see panels b, d vs. panels a, c in Fig. 7. Likewise, the imaginary part of the complex conjugate pair that arises at the merging point M 1 is small, and the oscillatory mode to which it gives birth will thus have a large period, as argued by Simonnet and Dijkstra (2002) for the oceanic gyre mode.
The same procedure (not shown) has been used for P 2 , traveling now back along the asymmetric branch, with α decreasing.We found again that the eigenvalue associated with the symmetry breaking merges with another eigenvalue to form an oscillatory mode.The onset of a jet mode at the pitchfork bifurcation P 2 is therefore quite similar to that at P 1 ; the only difference is in the direction of changing α: decreasing at P 2 vs. increasing at P 1 .In appendix A, we compute the atmospheric equilibria for other σ -values and we highlight the presence of other bifurcations.

Physical interpretation
It is possible to explain the Simonnet and Dijkstra (2002) scenario of LFV generation -with its succession of pitchfork → merging → Hopf bifurcation -in more familiar, physical terms: As α increases, the jet that lies close to the domain's north-south symmetry axis intensifies.At the first pitchfork bifurcation P 1 , the jet loses its stability and is subject the barotropic instability.Indeed, computing the Rayleigh criterion shows that β − u yy changes sign between 0 and P 1 , which is a necessary condition for barotropic instability to occur.mode.
The same procedure (not shown) has been used for P2, traveling now back along the asymmetric branch, with α decreasing.We found again that the eigenvalue associated with the symmetry breaking merges with another eigenvalue to form an oscillatory mode.The onset of a jet mode at the pitchfork bifurcation P2 is therefore quite similar to that at P1; the only difference is in the direction 340 of changing α: decreasing at P2 vs. increasing at P1.In appendix A, we compute the atmospheric equilibria for other σ-values and we highlight the presence of other bifurcations.

Physical interpretation
It is possible to explain the Simonnet and Dijkstra (2002) scenario of LFV generation -with its succession of pitchfork→merging→Hopf bifurcation -in more familiar, physical terms: As α 345 increases, the jet that lies close to the domain's north-south symmetry axis intensifies.At the first pitchfork bifurcation P1, the jet looses its stability and is subject the barotropic instability.Indeed, computing the Rayleigh criterion shows that β − uyy changes sign between 0 and P1, which is a necessary condition for barotropic instability to occur.
Increasing α further yields a PV surplus that can no longer be used to strengthen the jet since the Increasing α further yields a PV surplus that can no longer be used to strengthen the jet since the latter is now, in all likelihood, barotropically unstable.On the other hand, according to the results of the bifurcation study, the jet can undulate and remain stable along the branches of the pitchfork bifurcation.The maximum degree of undulation is given by the basin size and PV forcing pattern.As α increases further, the jet cannot sustain its distortion anymore and has to undergo another instability, which corresponds to the Hopf bifurcation highlighted in the previous subsection.
In other words, the pitchfork → merging → Hopf scenario found here for the first time as a generation mechanism of atmospheric LFV is closely connected to classical barotropic instability.It pursues, however, the further development of this instability into a highly nonlinear, chaotic regime.

Position of the pitchfork bifurcations
The pitchfork bifurcations described so far -and plotted in Fig. 5 above and in Figs.11 and A3 below -are however special cases of the broader picture illustrated in the regime diagram of Fig. 10.We tracked the pitchfork bifurcations in the (σ, α) plane of this figure by using -within a limited range of parameter values -the continuation algorithm described in Salinger et al. (2002).The number of curves in Fig. 10 gives an idea on the number of multiple equilibria for a given value of σ and α; it is, however, only a lower bound on their total number, due to the presence of isolated branches that have been encountered in the bifurcation diagrams so far.
The first pitchfork bifurcation, as σ increases, occurs at σ = 0.58.The higher the value of σ -and hence the longitudinal extent of the SST front -the more multiple equilibria we find.This reduces the range of σ -values of interest to that for which the number of equilibria is still reasonable.The number of pitchfork bifurcations at fixed σ does not appear to Nonlin.Processes Geophys., 19, 479-499, 2012 of the bifurcation study, the jet can undulate and remain stable along the branches of the pitchfork bifurcation.The maximum degree of undulation is given by the basin size and PV forcing pattern.
As α increases further, the jet cannot sustain its distortion anymore and has to undergo another instability, which corresponds to the Hopf bifurcation highlighted in the previous subsection.

355
In other words, the pitchfork→merging→Hopf scenario found here for the first time as a generation mechanism of atmospheric LFV is closely connected to classical barotropic instability.It pursues, however, the further development of this instability into a highly nonlinear, chaotic regime.18

Position of the pitchfork bifurcations
The pitchfork bifurcations described so far -and plotted in Fig. 5 above and in Figs.11 and 17 below -are but special cases of the broader picture illustrated in the regime diagram of Fig. 10.We tracked the pitchfork bifurcations in the (σ,α) plane of this figure by using -within a limited range of parameter values -the continuation algorithm described in Salinger et al. (2002).The number of curves in Fig. 10 gives an idea on the number of multiple equilibria for a given value of σ and α; it is, however, only a lower bound on their total number, due to the presence of isolated branches that have been encountered in the bifurcation diagrams so far.vary monotonically or in some other simple way with α.By plotting the model's total energy as a function of σ and α (not shown), one can verify that the curves in Fig. 10 never intersect: each curve is associated with a different energy level.

Atmospheric equilibria over an asymmetric ocean
There is at least one atmospheric equilibrium associated with each asymmetric equilibrium of the oceanic circulation.For such an oceanic circulation, the model's SST front has several meanders near the western boundary of the basin (cf.Fig. 3b).This situation is familiar from the Gulf Stream's flow downstream of Cape Hatteras (see, for instance, Lee andCornillon, 1996 or Dijkstra andGhil, 2005).
When the wind forcing σ increases, the number and size of the meanders increase, too.The corresponding atmospheric equilibrium is roughly aligned above the meanders of the SST front (not shown).When increasing α, we get an increasing number of atmospheric equilibria.This increase is due to the larger number of saddle-node bifurcations obtained when increasing σ , while using larger, but fixed values of α (not shown).
5 Temporal evolution of the atmosphere

Forcing by an antisymmetric ocean
We carried out several model integrations for σ = 0.67 and for different α-values; the corresponding equilibria appear in the bifurcation diagram of Fig. 11.For all these integrations, we chose to maintain a steady-state ocean circulation equal to the (unstable) antisymmetric solution for this σ -value, in order to isolate the atmospheric response, as done for σ = 0.6 and σ = 0.65 in Sect. 4. As before, we use the asymmetry measure E at to characterize the atmospheric state.In Fig. 12, we plot the estimated probability density function (PDF) of E at for seven different values of α, using in each case a 50 000-day long simulation, after the initial transients have died out.None of these curves is exactly symmetric with respect to zero, but the asymmetry is merely a sampling problem: a different choice of initial state can lead to a PDF that is either right-or left-skewed -perfect symmetry would require even longer time series.We denote by PDF α the estimated PDF associated with a given value of α.
At α = 0.07, there is no asymmetric equilibrium for the atmospheric flow (see again Fig. 11).The values of E at thus remain close to zero, as we can see from the curve PDF 0.07 in Fig. 12.A Fourier analysis of this time series (not shown) yields mainly two peaks, at periods of 7 and 14 days: they correspond to Rossby waves in the periodic β-channel.
For slightly larger values of α, the series takes on values that are strongly bimodal.The presence of multiple equilibria has a clear impact on the atmospheric model's time evolution: each peak of the PDF in Fig. 12 corresponds to a fixed point found in the bifurcation diagram of Fig. 11.We recover here the idea that the weather regimes are influenced by the existence of steady states, albeit unstable, of the evolution equation (Legras and Ghil, 1985;Ghil and Robertson, 2002).Indeed, the local maxima of the PDF are often interpreted as signatures of distinct weather regimes (Kondrashov et al., 2004).
As α increases, the range of the asymmetry E at in the atmospheric flows increases.These flows become also more irregular, and thus bimodality becomes less evident for α ≥ 0.12.Still, the gradually larger spread between the modes of PDF α as α increases in Fig. 12 is consistent with the fanning out of the asymmetric solution branches in Fig. 11  The SST forcing is similar to Fig. 3c.
local maxima of the PDF are often interpreted as signatures of distinct weather regimes (Kondrashov et al., 2004).
As α increases, the range of the asymmetry ∆Eat in the atmospheric flows increases.These 410 flows become also more irregular and thus bimodality becomes less evident for α ≥ 0.12.Still, the gradually larger spread between the modes of PDFα as α increases in Fig. 12 is consistent with the fanning out of the asymmetric solution branches in Fig. 11.

Forcing by an asymmetric ocean
We choose now to maintain a steady-state ocean circulation pattern that corresponds to a (likewise 415 unstable) asymmetric solution.When integrating the atmospheric model forced by the SST front of this oceanic solution, the climatological state of the atmospheric circulation is also shifted towards an asymmetric pattern, in which the average ∆Eat is no longer zero.This shift is explained by the fact that the asymmetric position of the SST front forces the atmospheric circulation in an asymmetric way, and does not inject the same positive and negative PV fluxes into the atmosphere.When the SST 420 front is deflected southward -which corresponds to ∆Eoc > 0, a strengthening of the subtropical 21 Fig. 11.Bifurcation diagram for the atmospheric equilibria at σ = 0.67; to be compared with the one obtained at σ = 0.65 and shown in Fig. 5c.The small crosses mark the position of saddlenode bifurcations; P 1 and P 3 are supercritical pitchfork bifurcations, while P 2 and P 4 are subcritical ones.The SST forcing is similar to Fig. 3c.

Forcing by an asymmetric ocean
We choose now to maintain a steady-state ocean circulation pattern that corresponds to a (likewise unstable) asymmetric solution.When integrating the atmospheric model forced by the SST front of this oceanic solution, the climatological state of the atmospheric circulation is also shifted towards an asymmetric pattern, in which the average E at is no longer zero.This shift is explained by the fact that the asymmetric position of the SST front forces the atmospheric circulation in an asymmetric way, and does not inject the same positive and negative PV fluxes into the atmosphere.When the SST front is deflected southward -which corresponds to E oc > 0, a strengthening of the subtropical gyre and a larger subpolar gyre -then the atmosphere's cyclonic circulation on the northern, cold side of the oceanic front is strengthened, while its anticyclonic circulation on the southern, warm side is weakened.The situation is reversed when the oceanic SST front is deflected northward.Indeed, the net PV input is positive when considering the mean PV flux over the entire basin.
The equivalent of Fig. 12, but for an asymmetric SST, is plotted in Fig. 13.The time series are obtained using the same set of parameters as in the previous set of runs but with the underlying SST shown in Fig. 3b.The PDFs are strikingly more peaked than those in Fig. 12.Indeed, this SST field maintains the atmosphere in an asymmetric state.The mean of the PDFs varies with the parameter α: as the forcing increases, the ocean imprints its asymmetric pattern on the atmosphere.
The climatological mean of a simulation of this type is different from the atmospheric equilibrium computed for the and for seven times series with the following α-values: 0.07, 0.08, 0.09, 0.10, 0.11, 0.12, and 0.13.The sharpest, most peaked PDF is for the smallest α-value, α = 0.07, the flattest is for the largest value, α = 0.13.The inserts are the mean states computed with (∆E at < 0; α = 0.1) on the left; (∆E at = 0; α = 0.07) for the peak of the PDF in the middle; and (∆E at > 0; α = 0.1) on the right.gyre and a larger subpolar gyre -then the atmosphere's cyclonic circulation on the northern, cold side of the oceanic front is strengthened, while its anticyclonic circulation on the southern, warm side is weakened.The situation is reversed when the oceanic SST front is deflected northward.Indeed, the net PV input is positive when considering the mean PV flux over the entire basin.

425
The equivalent of Fig. 12, but for an asymmetric SST, is plotted in Fig. 13.The time series are obtained using the same set of parameters as in the previous set of runs but with the underlying SST shown in Fig. 3b.The PDFs are strikingly more peaked than those in Fig. 12.Indeed, this SST field maintains the atmosphere in an asymmetric state.The mean of the PDFs varies with the parameter α: as the forcing increases, the ocean imprints its asymmetric pattern on the atmosphere.

430
The climatological mean of a simulation of this type is different from the atmospheric equilibrium computed for the corresponding oceanic pattern.Indeed, the atmospheric jet's meanders do not match those of the oceanic jet.The difference in size and intensity between the cyclonic and the 22 Fig. 12.Estimated normalized probability density functions (PDFs) of E at obtained for σ = 0.67 and for seven times series with the following α-values: 0.07, 0.08, 0.09, 0.10, 0.11, 0.12, and 0.13.The sharpest, most peaked PDF is for the smallest α-value, α = 0.07; the flattest is for the largest value, α = 0.13.The inserts are the mean states computed with ( E at < 0; α = 0.1) on the left; (α = 0.07) for the peak of the PDF in the middle; and ( E at > 0; α = 0.1) on the right.corresponding oceanic pattern.Indeed, the atmospheric jet's meanders do not match those of the oceanic jet.The difference in size and intensity between the cyclonic and the anticyclonic features is larger than in the steady state computed using the continuation method (not shown).Linking this PDF with the atmospheric steady states in the asymmetric case is a nontrivial task and would provide material for another paper.
In Fig. 14, we plot the difference between the SST field of the asymmetric state for which E oc > 0 and that of the antisymmetric state obtained for the same σ -value.This difference presents a tripole pattern that closely resembles the gyre-mode patterns, as shown in Fig. 4. The similarity results from the fact that the difference field in Fig. 14 is just a nonlinearly saturated version of the instability mode at the origin of the successive pitchfork bifurcations in Sect.3. Furthermore, the SST fields, whose difference is shown here, are essentially imposed by the corresponding stream function fields of the antisymmetric and asymmetric solutions.For this sign of the asymmetry, we observe a strengthening of the atmospheric cyclonic circulation north of the jet.

Summary
We used an oceanic and an atmospheric model with quasigeostrophic (QG), equivalent-barotropic dynamics in both fluids to search for multiple equilibria and low-frequency oscillations of the atmospheric and oceanic circulation.We studied separately the effects of the atmosphere on the ocean, as well as those of the ocean on the atmosphere.Doing so allowed us to focus on the atmosphere's intrinsic lowfrequency variability (LFV) when subject to forcing by the SST field.The originality of our ocean-atmosphere model lies in the parametrization of the Ekman boundary layer in the model's atmospheric component (cf.Feliks et al., 2004, cited as FGS in the main text).
With this parametrization -and using a horizontal resolution of 30 km in both the atmosphere and the oceanwe found that an oceanic sea surface temperature (SST) front does have a substantial influence on the dynamics of the free troposphere over and beyond the entire basin.This crucial observation provides a new framework for studying anticyclonic features is larger than in the steady state computed using the continuation method (not shown).Linking this PDF with the atmospheric steady states in the asymmetric case is a nontrivial 435 task and would provide material for another paper.
In Fig. 14, we plot the difference between the SST field of the asymmetric state for which ∆E oc > 0 and that of the antisymmetric state obtained for the same σ-value.This difference presents a tripole pattern that closely resembles the gyre-mode patterns, as shown in Fig. 4. The similarity results from the fact that the difference field in Fig. 14 is just a nonlinearly saturated version of 440 the instability mode at the origin of the successive pitchfork bifurcations in section 3; furthermore, the SST fields whose difference is shown here are essentially imposed by the corresponding stream function fields of the antisymmetric and asymmetric solutions.For this sign of the asymmetry, we observe a strengthening of the atmospheric cyclonic circulation north of the jet.

Summary
We used an oceanic and an atmospheric model with quasi-geostrophic (QG), equivalent-barotropic dynamics in both fluids to search for multiple equilibria and low-frequency oscillations of the atmo-23  14: Difference between the SST field of an asymmetric state for which ∆Eoc > 0 and that of the antisymmetric state, both obtained for σ = 0.9.Contours every 2 K; to be compared with Fig. 4. spheric and oceanic circulation.We studied separately the effects of the atmosphere on the ocean, as well as those of the ocean on the atmosphere.Doing so allowed us to focus on the atmosphere's intrinsic low-frequency variability (LFV) when subject to forcing by the SST field.The originality of our ocean-atmosphere model lies in the parametrization of the Ekman boundary layer in the model's atmospheric component, cf.Feliks et al. (2004, cited as FGS in the main text).
With this parametrization -and using a horizontal resolution of 30 km in both the atmosphere and the ocean -we found that an oceanic sea surface temperature (SST) front does have a substantial influence on the dynamics of the free troposphere over and beyond the entire basin.This crucial observation provides a new framework for studying the effects of mid-latitude SST anomalies on the mean atmospheric circulation and on its intraseasonal variability (Kushnir et al., 2002).
As a first step along this promising road, we studied the influence of changes in the mid-latitude oceanic, wind-driven circulation on the atmosphere.In the present, highly idealized setting, this study focused on the effects of varying two model parameters: the strength σ of the climatological wind forcing on the ocean, and the strength α of the ocean's SST forcing on the atmosphere.The parameter α is related to the height of the MABL by Eq. ( 3), while σ is prescribed in Eq. ( 4).The Fig. 14.Difference between the SST field of an asymmetric state for which E oc > 0 and that of the antisymmetric state, both obtained for σ = 0.9.Contours every 2 K; to be compared with Fig. 4. the effects of mid-latitude SST anomalies on the mean atmospheric circulation and on its intraseasonal variability (Kushnir et al., 2002).
As a first step along this promising road, we studied the influence of changes in the mid-latitude oceanic, wind-driven circulation on the atmosphere.In the present, highly idealized setting, this study focused on the effects of varying two model parameters: the strength σ of the climatological wind forcing on the ocean, and the strength α of the ocean's SST forcing on the atmosphere.The parameter α is related to the height of the MABL by Eq. (3), while σ is prescribed in Eq. ( 4).The effects of varying each of these parameters were covered in Sect. 3 for σ , and in Sects.4 and 5 for α.
The asymmetric branches are characterized -as in all similar QG, perfectly symmetric models -by a pairwise mirror symmetry.On these branches, several oscillatory modes destabilize the recirculation gyres near the merger of the boundary currents.The periodicity of the unstable modes is about 200 days for the fastest and 4 yr for the slowest, within the parameter range explored herein.Among the slowest modes is the interannual gyre mode first described by Jiang et al. (1995), which clearly differs in physical origin and spatial pattern from a classical, neutrally stable Rossby wave.In agreement with the careful diagnosis of Simonnet and Dijkstra (2002), this mode (Fig. 4) arises from the merging of two real eigenvalues, followed after a small σ -interval by a Hopf bifurcation; it is the smallness of the imaginary part of the resulting complex-conjugate mode, relative to its real part, that leads to the gyre mode's small, interannual frequency.
The asymmetric equilibria exhibit an oceanic SST front (Fig. 3b) of variable strength and extent, both of which depend on the intensity σ of the wind forcing.Sections 4 and 5 focused on the influence of such an SST pattern on the atmosphere.
In Sect. 4 and Appendix A, we plotted the bifurcation diagrams of the atmospheric circulation as α varies, for several σ -values (Figs. 5,11 and A3), while all other model parameters are equal to their values in Table 1.
As in the ocean, an antisymmetric state that is dominated by a westerly jet is present for all parameter ranges explored herein, but becomes unstable at fairly low σ -and α-values (Fig. 10).The asymmetric states that coexist with it are still pairwise mirror-symmetric, as for the ocean, but do not always arise from it via a pitchfork bifurcation, which makes their identification more difficult.
Focusing on the branch that originates from the first pitchfork bifurcation, at P 1 in Fig. 5, we found that -like in the oceanic case above -low-frequency instabilities arise from the merging of two real eigenvalues (Fig. 8).As before, one of these two is the one whose change of sign is responsible for the pitchfork bifurcation.The subsequent Hopf bifurcation, H 1 in Fig. 5, gives rise to an oscillation of the atmospheric jet on either side of its equilibrium position (Fig. 7).Therefore, we called these oscillatory modes -which are otherwise analogous to the gyre modes of the oceanic, doublegyre problem -the jet modes.
These jet modes do not always become unstable in the parameter range explored but, if they do, have a periodicity of about 70 days.Their physical nature and spatial structure, however, appear to be quite robust in all the model configurations explored herein.This instability is one possible explanation for LFV in the atmosphere.
Another type of temporal variability in our model's atmospheric component was studied in Sect.5: it corresponds to irregular jumps between the neighborhood of one or the other branch of asymmetric solutions, within a mirror-symmetric pair.The bimodal PDFs associated with this behavior were plotted in Fig. 12.
These two mechanisms -the jet mode and the bimodality associated with the coexistence of two unstable steady states -may help explain the dual character, episodic and oscillatory of atmospheric LFV, as highlighted by Ghil and Robertson (2002).On the one hand, the jet mode is responsible for low-frequency oscillations about the jet's equilibrium position.On the other hand, the steady states are responsible for the multimodal PDF that is often used to statistically diagnose weather regimes (Michelangeli and Vautard, 1998).

Discussion
Based on the results summarized above, we conclude that there is substantial similarity in the organization of the midlatitudes' oceanic and atmospheric circulation into multiple equilibria and oscillatory modes.In both cases, there is a strong link between symmetry breaking of the simplest, idealized flow pattern and the presence of a low-frequency instability: the gyre mode in the ocean and the jet mode in the atmosphere.A possible explanation of this similarity involves, of course, the fact that we consider very similar equations, as generally accepted in the theory of large-scale rotating flows (Gill, 1982;Ghil and Childress, 1987;Pedlosky, 1987).Still, the boundary conditions and the forcing are quite different, as they are for much more detailed models of atmospheric and oceanic flows.
Among the possible extensions to this study, we would like to understand how the atmospheric equilibria and oscillatory modes evolve when varying the configuration of the oceanic basin, so as to introduce realistic coastlines, or when considering baroclinic models for both the atmosphere and ocean.How the results might be affected by varying other parameters in Table 1 -in particular the oceanic and atmospheric Reynolds numbers R oc and R at -is likewise of considerable interest.Previous studies mentioned the impact of the Reynolds number on the oceanic and atmospheric circulation, and this aspect must be carefully investigated (see, for instance, Berloff et al., 2007a,b).
Another natural extension of this study would be to modify the MABL parametrization.As mentioned in the Introduction, one might want to take into account the vertical mixing, in order to parametrize the boundary layer height as a function of the underlying SST field, rather than accounting only for SST effects through the vertical velocity at the top of the MABL.Whether such changes in MABL height would noticeably affect the atmosphere's low-frequency variability (LFV) is an open question.
In a broader perspective, there are still several steps to be taken before being able to make a definitive connection between the instabilities identified in this study and the observed LFV over the Northern Hemisphere's ocean basins.The key step is to proceed to fully coupled, oceanatmosphere studies along the lines of the present study and of its predecessors.
We wonder whether this might not be also the case for the atmospheric problem under study here, as the oceanic SST front intensifies.
We thus consider again the atmospheric bifurcation diagram obtained at σ = 0.67, a slightly higher value than those explored in Sect.4, 0.6 ≤ σ ≤ 0.65; this diagram was shown in Fig. 11.The SST front at this σ -value is not very different from that at the value of σ = 0.65, used in Figs.5-9 (see again Fig. 3).The bifurcation diagram in Fig. 11, though, is quite different from the one in Fig. 5.
Up to and including the range 0.07 < α < 0.15, the two atmospheric bifurcation diagrams are fairly similar.We note, however, that the Hopf bifurcations H 1 and H 1 have migrated towards the pitchfork bifurcation P 1 , and H 2 and H 2 even further towards P 2 .In addition, there are two new pitchfork bifurcations: P 3 and P 4 , at α = 0.24 and α = 0.33, respectively.Each of these is associated with an additional real eigenvalue that becomes positive upon crossing the corresponding α-value.
In the same way as in Fig. 8 of Sect.4, we track here in Fig. A1 the evolution of the eigenvalue associated with the pitchfork bifurcation at P 3 along an asymmetric branch.Once more, we see that this eigenvalue merges with another real eigenvalue to form a complex conjugate pair -and thus an oscillating mode.From the pitchfork bifurcation at P 3 to the first saddle-node bifurcation at S 1 , the eigenvalue is negative.Between the two saddle-node bifurcations S 1 and S 2 , the eigenvalue is positive and the equilibrium is unstable.Beyond S 2 , the eigenvalue is negative again, only to merge, after a small α-interval, with another eigenvalue at M 1 .Note that all these transitions occur in a very narrow range of parameter values, 0.224 < α < 0.242.Furthermore, the bifurcation diagram in Fig. 11 shows that -unlike on the previous asymmetric branch that originates at P 2 -this mode does not become unstable for the range of α-values under scrutiny.
On the first pair of new asymmetric branches -the ones that originate at P 3 -we note the presence of two saddlenode bifurcations: S 1 and S 2 , and S 1 and S 2 , respectively.The spatial patterns of the eigenvectors associated with P 3 are plotted in Fig. A2.On none of the new branches in Fig. 11, though, does there appear to be a low-frequency Hopf bifurcation.
The eigenvector associated with the pitchfork bifurcation at P 3 appears in Fig. A2a, while those associated with the saddle-node bifurcations S 1 and S 2 are shown in Fig. A2b  and c, and those just before the merging point M 1 in Fig. A2d.In all four panels of Fig A2, the symmetric component of the eigenvector is dominant.The strong resemblance between Fig. A2b and c tends to prove that all pitchfork bifurcations in the atmospheric model arise due to the same eigenvector.
The eigenvector in Fig. A2a here, however, is less concentrated near the jet axis, and its stream function anomalies remain slightly larger even far from the jet axis.At M 1 , though, the eigenvector differs from the flow pattern obtained   The eigenvector in Fig. 16a here, however, is less concentrated near the jet axis, and its stream function anomalies remain slightly larger even far from the jet axis.At M1, though, the eigenvector differs from the flow pattern obtained in Fig 9c : wavenumber 1, in particular, is much stronger than in the previous case.Still, the largest amplitudes of the stream function are concentrated near the jet 580 axis and maintain an oscillation of the latter around its equilibrium position.
Just after the merging, and like in Fig. 8, it is again very difficult to follow the complex conjugate pair for larger values of α.At the merging point M1, the magnitude of the imaginary part is still very small with respect to the real part, although no low-frequency mode is generated.Moreover, the merging does not alter the spatial pattern of the eigenvector's real part.in Fig. 9c: wavenumber 1, in particular, is much stronger than in the previous case.Still, the largest amplitudes of the stream function are concentrated near the jet axis and maintain an oscillation of the latter around its equilibrium position.
Just after the merging, and like in Fig. 8, it is again very difficult to follow the complex conjugate pair for larger values of α.At the merging point M 1 , the magnitude of the imaginary part is still very small with respect to the real part, although no low-frequency mode is generated.Moreover, the merging does not alter the spatial pattern of the eigenvector's real part.

A2 Codimension-2 bifurcations
At σ 0.68, the connectivity of the asymmetric branches changes abruptly, as shown in Fig. A3: the branch that links the two pitchfork bifurcations, P 1 and P 2 , is severed, while the previously isolated branches in Figs. 5 and 11 become connected to the main antisymmetric branch.This type of sharp transition in branch connectivity is the signature of a Bogdanov-Takens (BT) bifurcation (Guckenheimer and Holmes, 1983).
For σ > 0.68, there is a range of σ -values for which we lost the two local equilibria (see also the description of the BT bifurcation in Kuznetsov (2004)).Several other features associated with the normal form of the BT bifurcation are present.These features include the merging and the splitting of two pairs of complex conjugate eigenvalues; the points where this occurs are marked by diamonds in Fig. A3b.
The gap between the two now separated branches increases when considering more intense SST fronts (not shown).We also conclude that, on this new branch, there

A2 Codimension-2 bifurcations
At σ ≃ 0.68, the connectivity of the asymmetric branches changes abruptly, as shown in Fig. 17: The branch that links the two pitchfork bifurcations, P 1 and P 2 , is severed, while the previously isolated branches in Figs. 5 and 11 become connected to the main antisymmetric branch.This type of sharp transition in branch connectivity is the signature of a Bogdanov-Takens (BT) bifurcation 590 (Guckenheimer and Holmes, 1983).
For σ > 0.68, there is a range of σ-values for which we lost the two local equilibria; see also the description of the BT bifurcation in Kuznetsov (2004).Several other features associated with the normal form of the BT bifurcation are present.These features include the merging and the splitting 29    exists an instability that is very similar to the jet mode patterns that were described in Sect. 4. On the branch emanating from P 3 , the jet modes appear at σ = 0.70, and become unstable as σ increases further (not shown).At σ = 0.70, this oscillatory mode has a period of 51 days.Its spatial pattern (not shown) is very similar to those already plotted in Fig. 7.

Appendix B
The marine atmospheric boundary layer (MABL) We use an analytical Ekman layer identical to that of FGS.The main approximation involved in their parametrization is to assume that the temperature in this layer is equal to the SST; this assumption corresponds to neutral conditions.The potential temperature θ in this layer is given by where T is the temperature, P the pressure, c v and c p the heat capacity of the air at constant volume and constant pressure, respectively, while the subscript ( where θ (x, y) = T oc (x, y) is the potential temperature anomaly in the MABL, which is equal to the SST anomaly at the same latitude and longitude.In the following, we omit the prime, since we concentrate on the anomalies only.
The zonal and meridional winds, u and v, in the MABL solve the equations: with f 0 the Coriolis parameter at 45 • N and k 0 the eddy viscosity.This coefficient is connected to the height of the layer by H E = π(2k 0 /f 0 ) 1/2 .By substituting the pressure from Eq. (B5) into these two equations, we can integrate them in z.
The only term that differs from the classical Ekman formulation is the temperature term in Eq. (B5) that expresses the thermodynamic pressure variation.The boundary conditions are that the horizontal velocity (u, v) vanishes at the ground and equals the geostrophic velocity (u G , v G ) at the top of the boundary layer: Using the continuity conditions at the top of the Ekman layer, we can integrate the Eqs.(B6) in the vertical so as to obtain an analytical expression of horizontal velocity in the boundary layer.Then, using the incompressibility equation that is integrated between the bottom and the top of the MABL, we obtain the vertical velocity w at its top, z = H E : The vertical velocity w(H E ) has two main components: one is the curl of the geostrophic wind speed and it is the traditional term associated with dissipation in the Ekman layer.The second term is proportional to the Laplacian of the temperature field in the MABL, and the coefficient of proportionality varies with the square of the height of the boundary layer, H E 2 .Thus, the choice of the value of H E is particularly important in determining the atmospheric dynamics and so we shall check this assumption in the main part of our study.Minobe et al. (2010) estimate the upward velocity above the Gulf Stream to be of the order of 0.05 Pa s −1 in winter and 0.01 Pa s −1 in summer.This order of magnitude is confirmed by Brachet et al. (2012).The correlation between the Laplacian of the SST field and the wind convergence has also been mentioned by Frankignoul et al. (2011) and by Shimada and Minobe (2011).According to Takatama et al. (2012), this mechanism explains, almost by itself, the entire vertical wind structure above the Gulf Stream.
atic model diagram.(a) Model cross-section; here ψ at is the atmospheric stream s the oceanic stream function, and T oc is the sea surface temperature (SST) field.The potential vorticity (PV) flux; see text for the definition of the other parameters.(b) e model, in which the solid box represents the contours of the oceanic basin, while delimits the contours of the atmospheric β-channel; this channel extends to the east oceanic basin.o atmospheric layers, of total height H at , with an Ekman layer of height H E and the re of height H at − H E .In this figure, we plot the fluxes between the three layers and urce terms of these fluxes.The term that forces the atmospheric vorticity is composed classical dissipation term, γ∇ 2 ψ at , and a forcing term proportional to the Laplacian d, α∇ 2 T oc .There is no external forcing of the atmospheric layer, such as a differential r heating.theocean is forced by a time-constant wind stress term σF, which mimics the clirface winds at midlatitudes, and subject to bottom friction with a constant relaxation dition, the SST is relaxed to an equilibrium profile chosen to resemble the temperature tropical oceanic basins; this SST forcing term is not shown in the figure.The z-axis , while the xand y-axes point east-and northward, respectively.

Fig. 1 .
Fig. 1.Schematic model diagram.(a) Model cross-section; here ψ at is the atmospheric stream function, ψ oc is the oceanic stream function, and T oc is the sea surface temperature (SST) field.The arrows mark a potential vorticity (PV) flux (see text for the definition of the other parameters).(b) Top view of the model, in which the solid box represents the contours of the oceanic basin, while the dashed box delimits the contours of the atmospheric β-channel; this channel extends to the east and west of the oceanic basin.
Fig. 3: (a,b) An asymmetric steady-state solution just below the first Hopf bifurcation H ′ 1 , with σ = 0.78.(a) Stream function ψ oc (x,y), with contours every 2 × 10 −2 ; stream function contours in nondimensional units, here and in similar figures.(b) SST field; contours every 2 K. (c) Antisymmetric SST field for the same value of σ.Here and in the following, positive contours are solid, negative contours are dashed.

Fig. 5 :
Fig.5: Bifurcation diagram for the atmospheric equilibria at σ = 0.65 .The asymmetry measure ∆Eat is computed in the same manner as in Fig.2for the ocean but as a function of the parameter α.All branches, whether stable or not, are shown as solid lines, for legibility.On the steady-state branches, we only plot -as open circles -the Hopf bifurcations for which the period exceeds 50 days.P1 is a supercritical pitchfork bifurcation and P2 is a subcritical one.The shape of the SST forcing used here is similar to Fig.3c.

Fig. 7 .
Fig. 7. Atmospheric stream function patterns spanning the oscillatory instability at the first two Hopf bifurcations.These patterns correspond to the (a, c) real and (b, d) imaginary part, respectively, of the atmospheric component of the eigenvectors associated with (a, b) H 1 and (c, d) H 2 ; H 1 and H 2 are marked with a circle in Fig. 5. Contours as in Fig. 6.

Fig. 8 :
Fig. 8: Leading real eigenvalue λ0 between the pitchfork bifurcation P1 and the merging point M1; point 2 is an intermediate point referred to in Fig. 9b.All three points are marked by solid triangles.

350Fig. 8 .
Fig. 8.Leading real eigenvalue λ 0 between the pitchfork bifurcation P 1 and the merging point M 1 ; point 2 is an intermediate point referred to in Fig. 9b.All three points are marked by solid triangles.

Fig. 9 :
Fig. 9: Real eigenvectors corresponding to the points identified in Fig. 8: (a) pitchfork bifurcation P 1 ; (b) intermediate point 2; and (c) merging of two real eigenvalues into a complex conjugate pair at M 1 .Contours as in Fig. 6.

Fig. 9 .
Fig. 9. Real eigenvectors corresponding to the points identified in Fig. 8: (a) pitchfork bifurcation P 1 ; (b) intermediate point 2; and (c) merging of two real eigenvalues into a complex conjugate pair at M 1 .Contours as in Fig. 6.

Fig. 10 :
Fig.10: Position of the pitchfork bifurcations in the (σ,α) plane.The 'S'-shaped curve is associated with the pitchfork bifurcations described in the previous subsections.The first pitchfork bifurcation, as σ increases, occurs at σ = 0.58.The higher the value of σ -and hence the longitudinal extent of the SST front -the more multiple equilibria we find.This reduces the range of σ-values of interest to that for which the number of equilibria is still reasonable.The number of pitchfork bifurcations at fixed σ does not appear to vary monotonically or in some other simple way with α.By plotting the model's total energy as a function of σ and α (not shown), one can verify that the curves in Fig.10never intersect: each curve is associated with a different energy level.

Fig. 10 .
Fig.10.Position of the pitchfork bifurcations in the (σ, α) plane.The S-shaped curve is associated with the pitchfork bifurcations described in the previous subsections.

Fig. 11 :
Fig. 11: Bifurcation diagram for the atmospheric equilibria at σ = 0.67; to be compared with the one obtained at σ = 0.65 and shown in Fig. 5.The small crosses mark the position of saddle-node bifurcations; P1 and P3 are supercritical pitchfork bifurcations, while P2 and P4 are subcritical ones.

Fig. 13 :
Fig. 13: Same as Fig. 12 but with an asymmetric SST field, shown in Fig. 3b.The insert is the mean state of the time series for α = 0.1.

Fig. 13 .
Fig. 13.Same as Fig. 12 but with an asymmetric SST field, shown in Fig. 3b.The insert is the mean state of the time series for α = 0.1.

Fig. 15 :
Fig. 15: Evolution of the leading real eigenvalue λ0 along the branch that arises from the pitchfork bifurcation P3 in Fig. 11 and on to the merging point M1; same conventions as in Fig. 8.The higher density of the points along certain segments of the curve in this figure reflects the continuation algorithm's slower convergence near those α-values.
Fig. A1.Evolution of the leading real eigenvalue λ 0 along the branch that arises from the pitchfork bifurcation P 3 in Fig.11and on to the merging point M 1 ; same conventions as in Fig.8.The higher density of the points along certain segments of the curve in this figure reflects the continuation algorithm's slower convergence near those α-values.

Fig. 17 :
Fig. 17: Bifurcation diagram for σ = 0.68: panel (b) is a blow-up of panel (a).The filled diamonds represent a splitting of a complex conjugate pair into two real eigenvalues, whereas the open diamonds mark the opposite situation.

Fig. A3 .
Fig. A3.Bifurcation diagram for σ = 0.68: (b) is a blow-up of (a).The filled diamonds represent a splitting of a complex conjugate pair into two real eigenvalues, whereas the open diamonds mark the opposite situation.
•) 0 indicates the reference state.Using the equation of state of an ideal gas, one hasln θ = − ln ρ + (1 − κ) ln P − ln R + κ ln P 0 ,(B2)with R the gas constant for the air and ρ its density.When differentiating this equation in the Boussinesq approximation(•) indicating a perturbed state with respect to the reference state (•) 0 .Using the Boussinesq approximation again, the hydrostatic equation is given hydrostatic equation over the height H E of the boundary layer, we find the expression of the pressure anomaly P at the height z: P (z) = P (H E ) + gρ 0 (z − H E )

Table 1 .
Model parameter values used in the standard case.

B. Deremble et al.: Multiple equilibria and oscillatory modes in a mid-latitude atmosphere model
Brachet et al. (2012)ble to the observed Gulf Stream or Kuroshio front -cf.FGS,Feliks et al. (2007),Feliks et al. (2011)andBrachet et al. (2012)-although the impact of eddies does affect this value in observed western boundary currents.