Multistability in a Coupled Ocean-Atmosphere Reduced Order Model: Non-linear Temperature Equations

Multistabilities were found in the ocean-atmosphere flow, in a reduced order ocean-atmosphere coupled model, when the non-linear temperature equations were solved numerically. In this paper we explain how the full non-linear Stefan-Bolzmann law was numerically implemented, and the resulting change to the system dynamics compared to the original model where these terms were linearised. Multiple stable solutions were found that display distinct ocean-atmosphere flows, as well as different Lyapunov stability properties. In addition, distinct Low Frequency Variability (LFV) behaviour was observed in stable attractors. We investigated the impact on these solutions of changing the magnitude of the ocean-atmospheric coupling, as well as the atmospheric emissivity to simulate an increasing green-house effect. Where multistabilities exist for fixed parameters, the possibility for tipping between solutions was investigated, but tipping did not occur in this version of the model where there is a constant solar forcing. This study was undertaken using a reduced-order quasi-geostrophic ocean-atmosphere model, consisting of two atmosphere layers, and one ocean layer, implemented in the Python programming language.

One source of non-linearity in the model equations comes from the long wave radiation emitted from the ocean and atmosphere, and modelled using the Stefan-Boltzmann law (σ B T 4 , where σ B is the Stefan Bolzmann constant).
In the MAOOAM model, the quartic radiation terms are linearised to simplify the projection of the equations onto a truncated Fourier expansion.This linearisation is justified by the fact that the perturbations in temperature are small in relation to the climatological reference temperatures (De Cruz et al., 2016).However, this linearisation removes the possibilities of non-linear interactions from the temperature terms.For this reason, this study investigates the impact on the LFV in the model, and the potential for multi-stabilities or bifurcations, when the temperature equations are not linearised.
We will show that removing the linearisation leads to multiple stable flow patterns in the atmosphere and ocean, for certain levels of ocean-atmosphere coupling and atmosphere emissivity.These flow patterns are qualitatively distinct and result in multiple average temperatures, for the same model parameters.They also present different cycle lengths and different dominant modes.
Section 2 describes the reduced order model used in this study, and Section 2.3 describes the modifications made to remove the requirement of linearising the temperature equations.This section also gives a description of the model configurations used in this study.The results are split into two sections, where we first look at the results from altering the ocean-atmosphere coupling (Section 3.1), and then the impact of atmospheric emissivity (Section 3.2).In each of the two results sections, we look at the impact of altering the given parameters on the stability and predictability of the system.Section 4 summarises the main results and discusses the general implications of these findings.

| Model Description
The qgs model (Demaeyer et al., 2020) is a reduced-order midlatitude climate model, with many different model configurations available.In the present work, we use the ocean-atmosphere model version where the atmospheric flow is obtained from a two-layer quasi-geostrophic flow defined on a β plane (Reinhold and Pierrehumbert, 1982).Similarly, the ocean streamfunctions are modelled using a quasi-geostrophic shallow-water model with a rigid lid (Pierini, 2011).The thermodynamic equation for the atmosphere and ocean temperatures are derived using an energy balance scheme proposed by Barsugli and Battisti (1998).The coupled ocean-atmosphere scheme used here was first introduced by Vannitsem et al. (Vannitsem et al., 2015).The atmosphercic variables are coupled through wind stress to the oceanic ones, driving the ocean circulation, which transports heat in the ocean.The ocean transfers heat with the atmosphere through radiative and direct heat coupling, which in turn impacts the atmospheric flow.
In this study we imposed a closed ocean basis (no-flux boundary conditions on all boundaries), and a channel atmosphere (no-flux boundary conditions on the north and south boundaries, and periodic boundary conditions at the west and east).We describe how these boundary conditions are implemented in Section 2.2.This version of the model with a closed ocean basin coupled to an atmosphere is called the Modular Arbitrary Order Ocean Atmosphere Model (MAOOAM) (De Cruz et al., 2016).
The governing partial differential equations (PDEs) for the atmosphere barotropic ψ a and baroclinic θ a streamfunctions and ocean streamfunctions ψ o are given as: where ψ a and ψ o are the atmosphere and ocean barotropic streamfunctions, and θ a are the atmosphere baroclinic streamfunctions.Vertical velocities are given by ω.
The ocean and atmosphere temperatures are derived from an energy balance model: where T a and T o are the atmosphere and ocean temperatures, γ a and γ o are the heat capacities of the atmosphere and ocean, σ is the static stability of the atmosphere (assumed constant), σ B is the Stefan-Boltzmann constant, R a and R o are the incoming solar radiation absorbed by the atmosphere and ocean, and ε is the atmospheric emissivity.
To reduce the number of variables, the atmosphere temperature variable T a is related to the baroclinic streamfunctions using the hydrostatic balance in pressure coordinates and the ideal gas law, providing the relationship

| Numerical Solution
The differential equations are projected onto basis modes, a procedure also known as Galerkin expansion.The basis modes are chosen to ensure that the boundary conditions described in the previous section are satisfied.This is done by stipulating that φ i (x, y ) = 0 for points (x, y ) on the boundary, and Sensible and turbulent heat exchange between ocean and atmosphere (Wm −2 K −1 ) R 287.058 Gas constant in dry air (Jkg −1 K −1 ) The key model parameters used in this study.Here C is a variable that we alter to investigate the impact of the ocean-atmosphere coupling.As in Charney and Straus (1980) we assume k = k d and as in Vannitsem (2015) that gravity g = 10ms −2 .
Three basis modes are of particular interest as they have real world analogies: • F 1 (x, y ) = √ 2 cos(y ) represents the solar insolation imbalance between the north and south.
• φ 2 (x, y ) = 2 sin(x /2) sin( 2y) is the double gyre, orientated so the peak is either to the north or south of the tough.
This loosely approximates the NAO, which is defined by the difference in surface pressure anomalies between northern and southern locations (often the Azores and Iceland) (Hurrell et al., 2003).The prevailing clockwise winds around the Azores high, and the counter clockwise winds around the northern low pressure can be broadly simulated by projecting the atmospheric streamfunction anomalies on this mode, thus simulating the impact on the wind and heat transport caused by the NAO.
The model variables are expanded using the basis modes.In previous studies using such energy balance models the temperature variables in the model are linearised around a fixed in time equilibria temperature: T (t , x, y ) = T 0 + δT (t , x, y ) to remove the quartic terms σ B T 4 (Vannitsem et al., 2015).This resulted in the temperatures being expressed as: δT o,i (t )φ i (x, y ). (7) The PDEs introduced in Equations (4, 5) are then projected onto these basis modes, using the inner product: This leads to 20 ordinary differential equations (ODEs) for the atmospheric streamfunctions, 10 for the barotropic and 10 for the baroclinic streamfunctions.In addition there are 16 ODEs in the ocean, 8 for the barotropic streamfunctions and 8 for the temperature anomaly.This leads to a total of 36 ODEs describing the model.

| Non-linear radiation terms
This study focuses on the change in the system dynamics caused by not linearising the radiation terms in the temperature equations of the MAOOAM model.This requires the reference temperature T a,0 and T o,0 to be time-varying quantities in the expansions shown in Equation 7. Therefore, to allow the average temperature across the atmosphere and ocean to change dynamically with time, we introduced two new basis modes, corresponding to constant spatial modes: F 0 (x, y ) = 1 and φ 0 (x, y ) = 1.These modes were added to the list of basis modes that were introduced in Section 2.2: . . .
These additional basis modes now allow the expansions, shown in Equation 7, to be given as T o (t , x, y ) = 8 i =0 T o,i (t )φ (x, y ), and similarly for the atmosphere.
Only the temperature equations are projected onto these additional basis modes, as these are the variables linearised in the original MAOOAM model.This means that we still have 10 ODEs for the atmospheric barotropic streamfunctions, and 8 ODEs for the oceanic barotropic streamfunctions.We obtain an additional ODE for the atmospheric baroclinic streamfunction (as this variable replaces the atmospheric temperature variables) and for the ocean temperature.This increases the total number of ODEs describing the system to 38.
We now introduce the two modified versions of the model that are used in the current study: Dynamic Equilibria: this version of the model includes the same linearisation as the Linearised version of the model, but the equilibrium temperature is dependant on time: T (t , x, y ) = T 0 (t ) + δT (t , x, y ).We will refer to this version as DE.
Non-Linear T4: This version does not involve the linearisation of the Stefan-Boltzmann law terms and the equations are projected directly onto the basis modes.This retains the quartic radiation terms.We will refer to this version of the model as T4.
This results in three model versions: • Linear Model (LM) • Dynamic Equilibria (DE) • Non-Linear (T4) More information about how these modifications were made can be found in the model manual (Demaeyer et al., 2022).See also the Supporting Information.

| RESULTS
The results section is split into two main parts.We first describe the system dynamics when the ocean atmosphere coupling is altered, and second we alter the atmospheric emissivity.We focus on values of these parameters that result in multiple distinct ocean-atmosphere flows, which are described by distinct attractors.Multiple stable attractors, for given parameter values, provide the possibility for model solutions to transition from one attractor to another, provided appropriate forcing is imposed.Such transitions between attractors could represent tipping points or abrupt changes or switching in flows.We are particularly interested in flows that present LFV, generated by the coupling between ocean and atmosphere variables, due to the increased predictability of the system dynamics that these solutions provide.

| Ocean-Atmosphere Coupling
This section presents the results where the magnitude of ocean-atmosphere coupling is varied for the three model versions.Increasing the ocean-atmosphere coupling has the effect of increasing the heat transfer between the ocean and atmosphere.This reduces the ocean temperature (which has a higher equilibria temperature than the atmosphere) and increases the atmospheric temperature, which has an impact on the baroclinic streamfunctions.The ocean temperature anomalies cause uneven temperature exchanges in the atmosphere that require the heat energy to be transported by atmospheric winds.At the same time, increasing the ocean-atmosphere coupling increases the friction felt by the atmosphere from the surface wind stresses with the ocean.In turn, this causes the atmospheric wind to have a greater impact on the movement of the sea temperature anomalies.Together, this reduces the fast moving timescales of the atmosphere through coupling to the ocean, which has a slower timescale relative to the atmosphere.
We alter the ocean-atmosphere coupling C by altering the following parameters: the strength of the oceanatmosphere coupling d , the ocean-atmosphere friction k d , the internal atmosphere friction k d , and the direct heat transfer between the ocean and atmosphere λ (Vannitsem, 2015).These parameters are controlled using a single friction coefficient C , where the relationship between C and the above parameters is given in Table 1.In this study we focus on values of C that are deemed to be within a realistic range (C ∈ [0.008kgm −2 s −1 , 0.02kgm −2 s −1 ]) for the real world coupling of the ocean and atmosphere (Houghton, 1986;Nese and Dutton, 1993;Vannitsem, 2017).

| State Space Probing
The effect of the ocean-atmosphere coupling on the system dynamics was investigated by fixing the coupling value C and running many trajectories with random initial conditions.Once these trajectories had appeared to settle onto an attractor, the initial transient section of the trajectories was discarded and we continued to run the trajectories for long run times (1 × 10 7 model days).This was done to ensure that the attractors remain stable and trajectories remain on the attractor.For these experiments, the atmosphere emissivity was set to ε = 0.7 as the model default (De Cruz et al., 2016).This process was repeated for different values of C .
The average temperatures of the ocean and atmosphere were calculated by projecting the trajectories, that were embedded within an attractor, onto the basis modes to obtain the temperature profiles for each time step.We then took the average across the spatial domain to obtain a single average temperature for each time step.Finally we took the average across time to obtain a single temperature, which represents the average temperature associated with a given attractor.This process is shown for the temperature of the ocean, with a similar process taken to calculate the average atmospheric temperatures: Where n o are the number of ocean modes, n x and n y are the number of spatial grid points being averaged across, and n τ is the number of time steps in the numerical solution.

| Multistabilies in Temperature -C
The average spatial and temporal temperatures are presented on a Temperature-C diagram.These diagrams are not bifurcation diagrams, as we only have information about the stable branches that we could find using the described method.The resulting figures for the average atmosphere and ocean temperatures, where the value of the ocean-atmosphere coupling is altered, are shown in Figure 1.In these figures, the branches are colour coded depending on the qualitative behaviour of the attractor.The figures show there are two intervals of C where there are multistabilies in the temperature.These intervals are approximately . These figures also show that there exists a single stable branch that bridges the two intervals (shown in blue).We plot the results of the dynamic equilibria (DE) and non-linear (T4) model runs on the same plot.
We found two differences between the T4 and DE runs.The first is that the difference in average temperatures between the stable branches in the DE runs are smaller than in the T4 runs.The small temperature differences between the attractors in the DE runs occurs due to the zeroth order temperature equations controlling the equilibria temperature (T a,0 , T o,0 ) having only one real stable solution.This means that the difference in average temperatures in each stable equilibrium is caused by only higher order terms.In the T4 model, there are additional terms in the zeroth order temperature equations that result in larger differences between the zeroth order temperatures.This difference in the equations occurs in the non-linear longwave radiation terms, and is sketched below for the ocean temperature equation.In the below expressions we introduce the tensor v i ,j ,k ,l ,m = F i , φ j , φ k , φ l , φ m .

Full Term
Zero Order Term The second difference between the model runs is that the T4 runs show a wider region of multistability in the . This is again assumed to be a result of higher order terms interacting in the T4 model equations, when compared to the DE runs, where the linearisation removes the higher order terms.
To investigate the properties of the attractors in the regions of multiple stability, the attractors were projected does not show LFV with any of the modes and this signifies that flows associated with this attractor do not present the same LFV or the same ocean dynamics as the other two attractors.The blue attractor becomes unstable for values C > 0.01275kgm −2 s −1 .This second bifurcation was found in Vannitsem (2017) using the linear version of the model, however the region of multistability in this interval was not found using the linear model.We have conducted long model runs to ensure the stability of both attractors in this interval and found that both attractors remain stable for at least 3 × 10 7 model days.This region of multistability was found in both the T4 and DE model runs, however it was found that the region of multistability is extended in the T4 model, when compared with the DE model.
Following the analysis of Vannitsem et al. (Vannitsem, 2017), we project the attractors onto the variables (ψ a,1 , ψ o,1 , T o,1 ) and (ψ a,1 , ψ o,2 , T o,2 ) to visualise the degree of ocean-atmosphere coupling in the attractors.In Figure 2 we show all three attractors, for the same values of C as before.In each image there is one attractor that stabilises around an unstable orbit that varies across all three variables over a decadal time scale.The image displaying (ψ a,1 , ψ o,2 , T o,2 ) shows that the attractor coloured in pink presents the same oscillating behaviour over the three variables as the attractor found in Vannitsem (2017).In addition, we also recover the attractor found by Vannitsem that does not present LFV (blue).The novel result here is that we have found an additional attractor (orange), for lower values of C , that also displays LFV in the coupled ocean-atmosphere system, but across the ocean-atmosphere modes (ψ a,1 , ψ o,1 , T o,1 ).
To visualise the resulting flow in the ocean, given the attractor behaviour, we have created videos showing the ocean streamfunction and temperature profiles, given the location in the projected 2D state-space, which can be found in the supplementary materials.The videos are named based on the colour coding used in the above plots.The

| Lyapunov Stability -C
To analyse the stability properties of the attractors that were found we calculated the Lyapunov properties of the stable attractors, focusing on values of C identified in the previous section where multiple stable attractors exist.For more details on calculating Lyapunov properties in coupled ocean-atmosphere models see Vannitsem and Lucarini (2016); Vannitsem (2017); Vannitsem et al. (2019); Vannitsem and Duan (2020).
We present the largest Lyapunov exponents (LLE) in Figure 3 (a), as a function of C , where again there is a multistability present in the two intervals identified.We have used the same colour coding as in Section 3.1.2to display which Lyapunov exponent is associated with which attractor.Approximately, the midlatitudes have a synoptic forecast time scale of less than two week (Lorenz, 1982), with larger scale process having a larger forecasting time (Lorenz, 1969).
At the synoptic scale, this would correspond to LLE of approximately 0.2-0.3days −1 .In the MAOOAM model there is a general decrease in the magnitude of the LLEs as C increases due to the increase in coupling between the ocean and atmosphere, resulting in the flow instabilities from the atmospheric dynamics being reduced (Vannitsem, 2017).However, it is clear that the orange branch displays significantly smaller LLEs than the other two branches, for relatively small values of C .Other studies have shown that for low values of C the magnitude of LLEs do decrease (Vannitsem, 2017), however this was observed for values C ≤ 0.0015kgm −2 s −1 , and we did not carry out runs for values of C this low due to these ocean-atmospheric coupling values being unrealistic in the real earth system.
In Figure 4  The extent to which coupling between the ocean and atmosphere is responsible for the LFV can be visualised using the variance of Covariant Lyapunov Vectors (CLVs) (Vannitsem and Lucarini, 2016).The CLVs provide a covariant basis of the tangent linear space of the system.In other words, the CLVs are vectors that form a basis and remain covariant with the flow, unlike forwards or backwards Lyapunov vectors (Kuptsov and Parlitz, 2012).Each covariant Lyapunov vector is stretched by the system dynamics by the corresponding local Lyapunov exponent.Each CLV is made up of 38 components, one for every variable of the system.The CLVs were calculated at each time step, and we took the variance of each of the 38 vector components across time, for each of the 38 vectors.The variance measures the variability of the CLV component in the direction of a given variable.Variables from the atmosphere and ocean that both have high variance for the same CLV index implies that these variables are interacting, or influencing each other.Therefore variables that are coupled with one another will present higher variance for the same CLV index (horizontal rows on the diagram).This allows us to visualise which variables are coupled and have a greater impact on the system dynamics.In Figure 5 we present a heatmap of the variance (log 10 scale) for the three distinct attractors.
These plots show the CLV index on the y-axis, and the model variables on the x-axis.The variables are ordered as: • Atmospheric barotropic streamfunctions ψ a (Index: 1-10) • Atmospheric baroclinic streamfunctions θ a (Index: 11-21) • Ocean barotropic streamfunctions ψ o (Index: 22-29) • Ocean temperature T o (Index: 30-38) As expected, the majority of the variance is seen in the atmosphere as these variables are the components that contribute to the fast timescale dynamics of the system.All three attractors present coupling for CLV indices 13-22, where the ocean temperature variables present similar variability to the atmosphere variables.The near zero Lyapunov exponents (CLV index 4-22) in general have a larger projection on the ocean variables (columns 22-38).All three attractors present similar dynamics for the indicies greater than 22, where the ocean components have low projection on the dynamics, implying that the large magnitude negative Lyapunov exponents are predominantly caused by the atmosphere dynamics.However, the attractors that present LFV have higher variance for the indices 1-12, implying Temperature-ε diagrams showing the averaged atmosphere (left) and ocean (right) temperatures.The different colours represent stable branches that present qualitatively different attractors.Here we only present the T4 model run results as the DE results produce the same results, and the temperature differences between the three branches of the DE runs are too small to be visible on these graphs.As in Figure 1, the attractors are projected onto the plane of ocean variables with the single gyre (left) and double gyre (right) for ε = 0.9.
that in these attractors the ocean temperature is having a stabilising impact on the atmosphere dynamics, and that there is a greater level of coupling between the ocean and atmosphere.In addition, the orange attractor (heat map on the left hand side) shows greater coupling between all four components of the model, as the ocean streamfunctions show higher variance for the indices 1-10.This explains the low magnitude of the positive Lyapunov exponents for this attractor as the unstable atmosphere components present high levels of coupling with the stable ocean.

| Emissivity
To simulate the impact of climate change on the ocean-atmosphere dynamics in the MAOOAM model, the emissivity ε is increased.Rising the emissivity acts as a proxy for rising levels of greenhouse gases, causing the atmosphere to 'trap' a larger proportion of the outgoing longwave radiation, thus increasing the average ocean and atmosphere temperatures.We picked a single value of ocean-atmosphere coupling C = 0.01175kgm −2 s −1 to investigate.In the previous section this value of C resulted in a single stable attractor, for ε = 0.7.

| Multistabilities in Temperatureε
By using Temperature-ε diagrams, shown in Figure 6, that present the average spatial and temporal temperatures of stable attractors found numerically, we can see that as the emissivity increases bifurcations occur at two values, leading to an additional two stable branches.We have coloured the stable branches to display the attractors that present qualitatively distinct behaviour.On these images we only present the results from the T4 model runs as each of the T4 and DE model runs presented the similar trajectory behaviour, but the resulting average temperature of the three distinct branches in the DE runs were too close to be visible on the graph.This is because of the higher order non-linear terms being removed in the linearised version, as explained in Section 3.1.2.
We have taken the three stable attractors that we found at ε = 0.9 and projected these onto the planes (ψ o,1 , T o,1 ) and (ψ o,2 , T o,2 ), and these are also shown in Figure 6.We can see that two of the attractors (shown in pink and blue) The three distinct attractors found at C = 0.01175kgm −2 s −1 and ε = 0.9 are projected onto (ψ a,1 , ψ o,1 , T o,1 ) and (ψ a,1 , ψ o,2 , T o,2 ), left and right respectively.The attractors are coloured to match the diagrams shown in Figure 6.
present the same behaviour as seen in Section 3.1, in addition, the highest temperature branch (shown in orange) qualitatively appears to have similarities with the attractor that showed a periodic behaviour with respect to (ψ o,1 , T o,1 ) that we also saw in the previous section.Therefore, increasing the value of ε appears to have the impact of providing stability to unstable branches.Figure 6 shows that the average temperatures in the atmosphere and ocean rise quickly as ε is increased, but in addition there are multistabilies that are separated by up to 1 • K for the same emissivity values.
As in Section 3.1, we visualise the level of ocean-atmosphere coupling by projecting the three stable attractors onto the first atmosphere mode, and the single and double gyre ocean modes.As with this previous section, we see that the orange attractor displays an oscillatory behaviour, coupling the barotropic streamfunctions and the single gyre variables, and the pink attractor displays LFV with the double gyre variables.Again, we see two distinct flows where there exists coupling between the ocean and atmosphere, over decadal time periods.These results show that in the MAOOAM model, rising emissivity leads to multistabilies in the ocean-atmosphere system, that were not present at low levels of emissivity.From our model runs we have not found examples of trajectories switching between the stable branches, however all of our model runs were undertaken using a constant solar forcing.To understand the robustness of the attractors to forcing, further model runs will have to be undertaken.
Similar to Section 3.1.2,we have produced videos to show the resulting ocean streamfunction and temperature behaviour given the attractors.These videos can be found in supplementary materials.The qualitative behaviour of each of the attractors identified in this section is similar to that of the corresponding attractors (those sharing the same colours) in Section 3.1.2.One minor difference in the results between the ocean-atmosphere coupling model runs and the emissivity model runs is that the LFV in the orange attractor over the first ocean mode φ o,1 is more clearly defined.This is shown by the projection of the attractor on the plane (ψ o,1 , T o,1 ), where the oscillating behaviour over these variables contains less noise and variation.This change in the orange attractor could be caused by the increase in ocean-atmosphere coupling between the two runs, where we used the value of C = 0.009kgm −2 s −1 in Section 3.1.2,and C = 0.01175kgm −2 s −1 in this section.

| Lyapunov Stabilityε
Following the format of Section 3.1.3,we present the LLEs in Figure 3 (b), where ε is varied, and we fix the oceanatmosphere coupling at C = 0.01175kgm −2 s −1 .As the value of ε is increased we see additional stable attractors appear, however the value of the LLEs on each stable branch does not alter with emissivity.This is because the emissivity has the impact of increasing the temperature of the layers evenly in space.The atmospheric layers in the model are driven from the meridional gradient in solar insolation, leading to baroclinic instability.In the current model set up the emissivity has no impact on this temperature gradient.To more closely model the expected outcomes of global heating, model runs should be undertaken where the rising emissivity reduces the temperature gradient between the Arctic and the equator (Francis and Vavrus, 2012;Rantanen et al., 2022).
We have presented the Lyapunov spectra of the three attractors at ε = 0.9, for C = 0.01175kgm −2 s −1 in Figure 4 (b).Interestingly, the orange attractor, with the lowest magnitude LLE, shows a significantly lower number of near zero Lyapunov exponents, compared to the other two attractors and also the orange attractor presented in Section 3.1.3.There is a clear drop in the magnitude of the Lyapunov exponents at the index 18.This is an interesting result as it implies that the Lyapunov dimension of this attractor is significantly lower than the other attractors.This difference between the Lyapunov spectra of the orange attractors in Figure 4 (a) and (b) could be caused by the increase in ocean-atmosphere coupling between the two model runs, from C = 0.009kgm −2 s −1 to C = 0.0175kgm −2 s −1 in this section.

| DISCUSSION
We have described the novel multistabilities found in a reduced order atmosphere-ocean model, resulting from not linearising the longwave radiation terms With one attractor producing the double gyre behaviour, similar to what is observed over the North Atlantic, and the other attractor displaying a more complex flow, where the main relationships are with the first ocean mode φ 1 and the fifth ocean mode φ 5 .This attractor displayed the greatest degree of coupling between the ocean and atmosphere, comparing with the other stable attractors, where all four variables are coupled.This leads to significantly lower positive Lyapunov exponents, which implies that this attractor would have a longer forecasting window.Further studies will need to be done to find if this attractor describes a real-world ocean-atmosphere flow.
A key reason for undertaking this study, and not linearising the longwave radiation terms, was to investigate the potential of tipping between multistabilities.While we found cases of distinct attractors for the same parameter values, we could not produce trajectories that switched intermittently between the stable branches, though it is not possible to rule out the possibility of such trajectories existing.However, in the current model setup all external forcings are stable with respect to time.To test the robustness of the stability of each attractor model runs could be undertaken, where stochastic forcing or perturbations are included, to see if there is the potential for noise induced tipping between the stable branches.Another potential source of tipping could come from periodic cycles, such as the annual solar cycle as implemented in a similar linearised model (Vannitsem, 2017).
With rising global temperatures, investigating the possibility and impact of tipping points in the climate is of great importance.Understanding how rising temperatures could impact existing multi-decadal patterns in the climate could lead to a better understanding of how established climate patterns may look in the future.Or if there is the possibility of abrupt transitions from one regime to another.This paper has introduced a modified model that produces both multistabilities, as well as attractors that present LFV, that become stable for rising emissivity.These properties are of interest as it facilitates the study of attractors that allow forecasting well beyond the atmospheric Lyapunov time, as well as the potential of tipping from one attractor to the other.These aspects will be investigated in the future.
x on the boundary and F i (x, y ) = 0 for y on the boundary.In this study we use 10 basis modes for the atmosphere and 8 for the ocean, as inVannitsem (2017).These are set on a domain of x ∈ [0, 2π ], y ∈ [0, π ].The atmosphere F i and ocean φ i modes are given below: × 10 −11The meridional gradient of the Coriolis parameter at a given latitude (m −1 s −1 ) Atmosphere-surface friction (s −1 ) k d gC /∆p Internal atmosphere friction (s −1 ) r 1 × 10 −7 Ocean bottom Rayleigh friction (s −1 ) Depth of the ocean layer (m) d C /(ρ o h) Coefficient of mechanical ocean-atmosphere coupling (s −1 ) γ a 1 × 10 7 Specific heat capacity of the atmosphere (Jm −2 K) γ o 5.46 × 10 8 Specific heat capacity of the ocean (Jm −2 K −1 ) σ 0.2 Static stability of the atmosphere λ 1004C

F
I G U R E 1 The Temperature-C diagrams for the averaged atmosphere (left), and ocean (right) temperatures for the DE model (crosses), and the T4 mode (dots), of the stable attractors.The branches were colour coded by qualitatively investigating the behaviour of each attractor.The attractors are projected onto the planes (ψ o,1 , T o,1 ) (left) and (ψ o,2 , T o,2 ) (right), to display how the behaviour differs between the three stable attractors.The projection of the orange attractor is shown for C = 0.009kgm −2 s −1 and the other two projections (blue and pink) of the attractors are shown for C = 0.0125kgm −2 s −1 .
The distinct attractors for the T4 runs are projected onto (ψ a,1 , ψ o,1 , T o,1 ) (left) and (ψ a,1 , ψ o,2 , T o,2 ) (right), to present the level of ocean-atmosphere coupling.The orange attractor is shown for C = 0.009kgm −2 s −1 and the other two attractors are shown for C = 0.0125kgm −2 s −1 .Here we only present the T4 runs as the dynamics over these three variables are identical between the T4 and DE model runs.
onto the ocean variables (ψ o,1 , T o,1 ) and (ψ o,2 , T o,2 ) to visualise the behaviour.Looking at these projections we can see that there are three qualitatively distinct attractor behaviours, which correspond to three distinct flow behaviours in the ocean.We project all three attractors onto the same figure for comparison, shown in Figure1, with the orange attractor shown for C = 0.009kgm −2 s −1 , and the other two attractors shown for C = 0.0125kgm −2 s −1 .In the interval C ∈ [0.008kgm −2 s −1 , 0.009kgm −2 s −1 ], there exists two attractors, where the attractor shown in orange displays LFV with respect to the variables (ψ o,1 , T o,1 ).The orange attractor becomes unstable when C > 0.009kgm −2 s −1 .This multistability has not been found in other studies using the linear MAOOAM model.The LFV in the first ocean modes suggests that the variability in the ocean temperature and streamfunctions is prominently impacted by the single gyre oscillation in this case.The pink attractor, which becomes stable for values C > 0.012kgm −2 s −1 , presents LFV over the variables (ψ o,2 , T o,2 ), (a) LLE of the DE and T4 model runs, for varying C and ε = 0.7.(b) LLE of the DE and T4 model runs, for C = 0.01175kgm −2 s −1 and varying ε.F I G U R E 3 The Largest Lyapunov Exponents (LLE), shown in days −1 , where the results are colour coded based on the attractor behaviour, as described in Sections 3.1.2and 3.2.1.signifying that the dynamics of the flow are impacted by the double gyre dynamics.The other attractor (blue) present attractor coloured in blue shows a persistent positive value for the T o,2 variable, leading to a persistent double gyre temperature anomaly in the ocean temperature.In addition no clear LFV is present in the videos.The pink attractor shows oscillating behaviour over the variables T o,2 and ψ o,2 with timescales of approximately 70 model years, leading to transitions between a double and quadruple gyre temperature anomaly in the ocean.Similarly, a clear oscillation is seen in the ocean streamfunction where the most prominent variables oscillate between ψ o,2 and ψ o,6 over the same time period.Lastly, in the case of the orange attractor, no clear double gyre appears in the ocean temperature profile, and the LFV instead manifests in the first ocean mode φ o,1 .This leads to LFV with a timescale of 80-100 model years over the temperature variables T o,1 , T o,5 , and T o,7 .(a) Lyapunov spectra where the orange attractor is shown for C = 0.009kgm −2 s −1 , and the other two for C = 0.0125kgm −2 s −1 , where ε = 0.7.(b) Lyapunov spectra for ε = 0.85 and C = 0.01175kgm −2 s −1 F I G U R E 4 The Lyapunov spectra, shown in days −1 , for the three distinct attractors found are plotted to compare the magnitude of the largest Lyapunov exponent λ 1 , and the number of near zero Lyapunov exponents.
(a)  we present the Lyapunov spectra of the three distinct attractors identified while varying the level of ocean-atmosphere coupling.The figure shows that the two attractors at C = 0.0125kgm −2 s −1 (blue and pink), both have 18 near zero Lyapunov exponents, 17 negative exponents, and three positive exponents.However the positive Lyapunov exponents are approximately double in one of the attractors compared with the other, showing that the rate of divergence of initial conditions will be much greater in the blue attractor.The novel attractor found in this study (orange), presents significantly smaller magnitude positive Lyapunov exponents than the other attractors, and has a lower number of near zero Lyapunov exponents, implying that the attractor exists on a lower dimensional chaotic manifold than the other attractors.The smaller positive magnitude Lyapunov exponents for two of the three attractors occurs due to these attractors existing around unstable periodic orbits that produce the LFV on decadal timescales and therefore increases the predictability.F I G U R E 5 The variance of the CLVs (shown on a log 10 scale), for each of the three distinct attractors shown in Figure 1.The attractors are designated by the colour of their title, which corresponds to the previous sections.The orange attractor CLVs are calculated for ε = 0.7, and C = 0.009kgm −2 s −1 , and the other two attractors' CLVs are calculated with ε = 0.7 and C = 0.0125kgm −2 s −1 .
. The modifications made to the MAOOAM model have resulted in several features that were not present in the original linearised version with fixed reference temperature.This study introduced two new versions of the model (the dynamic equilibria and non-linear versions), and compared the results of these versions with the existing linear model.The properties of the new attractor found, as well as the region of multistability, were analysed qualitatively and by using the Lyapunov properties of the attractors.We have demonstrated in the reduced order ocean-atmosphere model, MAOOAM, that by modifying the linearisation of the longwave radiation terms we can obtain three qualitatively distinct stable attractors, with intervals of multistability for certain parameters.Interestingly, the dynamic equilibria version of the model, which includes the same linearisation as the original model but allows the zeroth order equilibria temperature to change with time, presents similar dynamics as the fully non-linear version for most parameter values.All three distinct attractors can be obtained in the dynamic equilibria (DE) version, where there are two attractors that present LFV while representing largely different coupled ocean-atmosphere flows.In addition, the DE version of the model can be run in the same length of time as the fixed reference temperature version, which is almost an order of magnitude faster than the non-linear version of the model.Two of the distinct attractors present low frequency variability (LFV) behaviour.One of the attractors, which displays LFV over the second ocean mode φ 2 , has been identified in the linear model.In this study we found an additional attractor that displays LFV over the first ocean mode φ 1 , which has not been identified in the linear model.This attractor has a longer time scale (∼ 80 − 100 years) compared with the first attractor (which displays a timescale of approximately 70 years).In addition the two attractors display marked differences in the ocean and atmosphere flows.