Exploring the Lyapunov instability properties of highdimensional atmospheric and climate models
Abstract
The stability properties of intermediateorder climate models are investigated by computing their Lyapunov exponents (LEs). The two models considered are PUMA (Portable University Model of the Atmosphere), a primitiveequation simple general circulation model, and MAOOAM (Modular ArbitraryOrder OceanAtmosphere Model), a quasigeostrophic coupled oceanatmosphere model on a plane. We wish to to investigate the effect of the different levels of filtering on the instabilities and dynamics of the atmospheric flows. Moreover, we assess the impact of the oceanic coupling, the dissipation scheme and the resolution on the spectra of LEs.
The PUMA Lyapunov spectrum is computed for two different values of the meridional temperature gradient defining the Newtonian forcing to the temperature field. The increase of the gradient gives rise to a higher baroclinicity and stronger instabilities, corresponding to a larger dimension of the unstable manifold and a larger first LE. The KaplanYorke dimension of the attractor increases as well. The convergence rate of the rate functional for the large deviation law of the finitetime Lyapunov exponents (FTLEs) is fast for all exponents, which can be interpreted as resulting from the absence of a clearcut atmospheric timescale separation in such a model.
The MAOOAM spectra show that the dominant atmospheric instability is correctly represented even at low resolutions. However, the dynamics of the central manifold, which is mostly associated to the ocean dynamics, is not fully resolved because of its associated long time scales, even at intermediate orders. As expected, increasing the mechanical atmosphereocean coupling coefficient or introducing a turbulent diffusion parametrization reduces the KaplanYorke dimension and KolmogorovSinai entropy. In all considered configurations, it is possible to robustly define large deviations laws describing the statistics of the FTLEs corresponding to the strongly damped modes, while the opposite holds for nearzero LEs and, somewhat unexpectedly, also for the positive LEs.
This paper highlights the need to investigate the natural variability of the atmosphereocean coupled dynamics by associating rate of growth and decay of perturbations to the physical modes described using the formalism of the covariant Lyapunov vectors and to consider long integrations in order to disentangle the dynamical processes occurring at all time scales.
[1]LesleyDe Cruz \Author[2]SebastianSchubert \Author[1]JonathanDemaeyer \Author[2,3,4]ValerioLucarini \Author[1]StéphaneVannitsem
1]Royal Meteorological Institute of Belgium, Brussels, Belgium 2]Meteorological Institute, CEN, University Of Hamburg, Germany 3]Department of Mathematics and Statistics, University of Reading, UK 4]Centre for the Mathematics of Planet Earth, University of Reading, UK
Lesley De Cruz (lesley.decruz@meteo.be)
The dynamics of the atmosphere and the climate system is characterised by the property of sensitivity to initial states (Kalnay, 2003). This feature implies that any small errors in the initial conditions will progressively amplify until the forecast becomes useless, or in other words cannot be distinguished from any random state taken from the climatology of the system. This property was already recognised in the early developments of weather forecasts (Thompson, 1957) and was associated with the nonlinear nature of deterministic dynamical systems by Lorenz (1963). These pioneering works sowed the seeds for the development of predictability theories for the atmosphere and climate, and for important progress in the context of dynamical systems, in particular the development of chaos theory (Eckmann and Ruelle, 1985). This sensitivity property affects not only the dynamics of errors in the initial conditions but also the errors that are present either in the model parametrizations, known as model errors, or in the boundary conditions (Nicolis, 2007; Nicolis et al., 2009). Clarifying the nature of this sensitivity is therefore crucial in the perspective of improving forecasts at short, medium and long term (Vannitsem, 2017)
The property of sensitivity to initial conditions in deterministic dynamical systems is often evaluated by computing the Lyapunov exponents that correspond to the asymptotic rates of amplification or decay of infinitesimally small perturbations, e.g. Eckmann and Ruelle (1985); Ott (2002) and Cencini et al. (2010). A system is chaotic if it possesses at least one positive Lyapunov exponent. Since the eighties many dynamical systems in various domains of science have been analysed from this perspective. This has revealed the presence of chaos in systems ranging from the fields of chemistry and biology to turbulence, e.g. Yamada and Ohkitani (1988); Gallez and Babloyantz (1991); Manneville (1995) and Sprott (2010). In the early days the investigations essentially dealt with loworder systems, but later the scope was broadened to include spatially distributed systems with a large number of degrees of freedom, in coupled maps, e.g. Nicolis et al. (1992); Vannitsem and Nicolis (1996); Cencini et al. (2010) and Yang and Radons (2013), and in partial differential equations, e.g. Manneville (1985, 1995); Vannitsem and Nicolis (1994) and Yang and Radons (2013). Recently, Lyapunov analysis was the subject of a special issue edited by Cencini and Ginelli (2013).
In parallel to these investigations in the context of basic sciences, several attempts to compute Lyapunov exponents in the context of meteorological and climate models have been made, see Vannitsem (2017), in particular in intermediateorder atmospheric quasigeostrophic models (with O(1000) variables) (Legras and Ghil, 1985; Vannitsem and Nicolis, 1997; Lucarini et al., 2007; Schubert and Lucarini, 2015, 2016). These analyses indicate that if realistic boundary conditions and forcings are imposed on the model under investigation, the number of positive exponents is high, which implies that the solution for the atmosphere lives on a highdimensional attractor. This suggests at first sight that the number of degrees of freedom needed to describe the dynamics is high and cannot be reduced to a loworder system.
The atmosphere cannot be treated as an autonomous system, as it interacts with other components of the climate system. These other components are characterised by longer time scales of motions. They are typically less intensely affected by some of the physical processes responsible for atmospheric instabilities, and most notably convective and baroclinic instability. Moreover, the energetics of the atmosphere is mainly driven by thermodynamic processes that are dominated by the inhomogeneous absorption of solar radiation. The oceanic circulation, by contrast, is mostly mechanically driven by atmospheric winds (Lucarini et al., 2014). This raises the question as to the impact of the coupling to other subdomains of the climate system: are the other subdomains of the climate system stabilising the atmosphere or not? Vannitsem et al. (2015) partly addressed this question in the context of coupled loworder oceanatmosphere systems. They found that the presence of the ocean and its exchanges (heat and momentum) with the atmosphere can drastically reduce the instability properties of the flow, confirming earlier results of Nese and Dutton (1993). As discussed below, the role of the ocean in modulating and impacting atmospheric instabilities is far from trivial.
Yet the problem of the predictability (in terms of Lyapunov instability) of the fullscale climate system including the different subdomains is still open. Recently a new coupled oceanatmosphere model was developed that could help answer key questions on the predictability properties of this type of system (De Cruz et al., 2016). The model was coined MAOOAM for Modular ArbitraryOrder OceanAtmosphere Model. The modular design of MAOOAM allows one to easily explore different model parameters and resolutions. In particular, the coupling strength between the ocean and the atmosphere should modify the predictability properties of the flow as illustrated in (Vannitsem et al., 2015). Moreover, the model resolution is also expected to play an important role in the instability properties of the flow as discussed in (Lucarini et al., 2007) in the context of an atmospheric model.
0.1 The properties of the tangent space
As originally envisioned by Ruelle (1979), it is possible to associate to each Lyapunov exponent a corresponding infinitesimal perturbation that covaries with the orbit that grows or decays asymptotically with the rate given by the corresponding exponent. These physical modes are usually referred to as covariant Lyapunov vectors (CLVs). The application of such a formalism to explore the properties of the tangent space was pioneered by Legras and Vautard (1995) and Trevisan and Pancotti (1998), before Ginelli et al. (2007) and Wolfe and Samelson (2007) provided efficient algorithms to compute them for highdimensional systems. The CLVs have been used to study e.g. spatiotemporal chaos (Pazó et al., 2008, 2010), RayleighBénard convection (Xu and Paul, 2016), and the dynamics of the midlatitudes atmosphere in the quasigeostrophic (QG) approximation (Schubert and Lucarini, 2015, 2016). Schubert and Lucarini (2015, 2016) have also underlined that CLVs allow for generalising the classic normal mode instability of fixed stationary states to the case of chaotic background state (e.g. Charney, 1947; Eady, 1949; Pedlosky, 1964). Trevisan et al. (2010) and Carrassi et al. (2008) also showed that performing data assimilation on the unstable manifold spanned by the CLVs corresponding to positive Lyapunov exponents is extremely efficient because it allows one to focus on the portion of the tangent space supporting the growth of errors.
Additionally, CLVs allow for understanding the properties of the tangent space and assess the hyperbolicity of the system, through the analysis of the statistics of the angles between the stable and unstable tangent manifolds across the attractor. These angles should always be bounded away from 0 or in the ideal case of uniform hyperbolicity. This point of view complements the investigation of the statistical properties of FTLEs: the probability density functions of the FTLEs whose long term averages correspond to positive exponents do not cross zero in the case of uniform hyperbolicity. Note that the uniform hyperbolicity is key to defining the structural stability of a chaotic dynamical system and provides, through the chaotic hypothesis by Gallavotti and Cohen (1995), an important working hypothesis for constructing the statistical mechanics of highdimensional chaotic systems, even in the case that such a system is not uniformly hyperbolic. Note that uniform hyperbolicity also allows for establishing a rigorous response theory for chaotic dynamical systems (Ruelle, 2009), which has also been shown to apply well in complex systems where there is no reason to believe that such stringent condition on the tangent space is obeyed; see, e.g., Lucarini et al. (2017).
0.2 Multiscale Properties
As well known, geophysical fluid dynamical (GFD) systems are characterised by relevant processes on multiple spatial and temporal scales of motion (Schneider, 2006; Vallis, 2006). These scales of motion can be isolated by assuming dominant dynamical balances and performing corresponding asymptotic expansions of the dynamical equations (Klein, 2010). A possible way to look at the signature of such a diverse range of dynamical processes in a nonlinear, chaotic setting can be found by considering the general idea proposed by Gallavotti and Lucarini (2014) according to which one can expect to find that LEs corresponding to smaller time scales are associated to CLVs characterised by small spatial scales. By looking at the properties of the structure of each CLV, one should ideally be able to understand what kind of dynamical processes (e.g. QG vs. mesoscale) are mainly responsible for such a physical mode.
The problem becomes particularly interesting when considering the coupling of two subdomains with vastly different time scales as done in the case of a loworder coupled oceanatmosphere system in Vannitsem and Lucarini (2016). Three different manifolds were isolated in the model, the usual (highly) unstable manifold mainly associated with the dynamics of the atmosphere, a highly dissipative manifold also mainly associated with the dynamics of the atmosphere, and an extremely weakly (un) stable manifold, that will be here referred to as the central manifold, essentially dominated by the dynamics and thermodynamics of the ocean but coupled to the atmosphere as well. The presence of a nontrivial central manifold is typical of the socalled partially hyperbolic systems (Pesin, 2004). The CLVs corresponding to the central manifold are geometrically quasidegenerate, so that errors propagate easily between the various modes and impact both the atmosphere and the ocean. The corresponding FTLEs are strongly correlated and each have a rather slow decay of correlations, so that large deviation laws cannot be effectively estimated (Kifer, 1990; Touchette, 2009; Pazó et al., 2013; Laffargue et al., 2013). A particular consequence of this feature is that errors affecting the central manifold display a complex superexponential behaviour. The question is therefore what should be the resolution of the coupled atmosphereocean model and what should be the time of observation such that a better separation emerges between such modes.
0.3 This paper
In this paper we wish to provide some first steps of a wider research programme aimed at performing a systematic investigation of the properties of the tangent space of GFD systems in a turbulent regime of motion. A first objective is to gain a better understanding of the multiscale properties of the dynamics and of the energy exchanges occurring across such scales. Furthermore, this programme aims at understanding the relevance of violations to the uniform hyperbolicity conditions in terms of predictability on different time scales, including the response – in a statistical mechanical sense – of the system to static and timedependent perturbations.
In the present manuscript, we explore for the first time the Lyapunov spectra of a primitiveequation model, PUMA, and of the intermediateorder coupled oceanatmosphere system, MAOOAM. The first model is characterised by the presence of multiple scales of motions resulting from the fact that ageostrophic motions are not filtered, as opposed to the QG case (Klein, 2010). Instead, in the second model the multiscale properties come from the fact that the represented two geophysical fluids have largely different internal time scales.
For PUMA, we consider the first 200 Lyapunov exponents for two different meridional temperature gradients. We study the properties of the Lyapunov spectrum and on the estimates of the KaplanYorke dimension and KolmogorovSinai entropy (Eckmann and Ruelle, 1985). In the case of MAOOAM, we investigate the role of dissipation introduced in the model (linear friction and effective diffusion) and the impact of the resolution of the models. For both models, the existence of large deviation laws of the FTLEs is tested.
In Section 1, the two models are described and Section 2 is devoted to a brief description of the Lyapunov instability analysis and the experimental setups. Section 3 summarises the results obtained so far and in Section 3.2 we present our future programme, which aims to clarify the instabilities of highresolution systems.
1 Model description
1.1 Puma
The Portable University Model of the Atmosphere (PUMA) was introduced by Fraedrich et al. (1998). The intent of its developers was to design a model that gets close to stateoftheart circulation models and at the same time is still easy to use in teaching and research by single scientists. PUMA is the dynamical core of the Planet Simulator (PLASIM) which is a fully coupled climate model of intermediate complexity. PLASIM has been frequently used to study storm tracks (Fraedrich and Kirk, 2005), explore the sensitivity of the sea ice albedo bifurcation (Boschi et al., 2013) or create large ensembles for climate change experiments allowing for interesting applications in economic modelling (Holden et al., 2013) or to assess the feasibility of linear response theory (Ragone et al., 2015).
Let us briefly summarise the equations of motion of PUMA and how the model is integrated in time. For further details we refer the reader to the PUMA User’s Guide (Fraedrich et al., 1998). PUMA solves the primitive equations, which are derived from the NavierStokes equation on a rotating sphere by assuming approximate hydrostatic balance. This means that (convective) motions which are characterised by a vertical acceleration with a size comparable to gravity are filtered out (Holton, 2004). The prognostic equations as written in PUMA’s code have four prognostic fields, the relative vorticity , the divergence , the temperature and the logarithm of the surface pressure . These equations are
(1)  
(2)  
(3)  
(4) 
where the vorticity is defined as and the divergence is defined as . Additionally, one takes into account the hydrostatic relation
(5) 
and is defined as with . Some abbreviations have been used:
with . The variables used in this equation can be found in Table 1.
PUMA is forced by Newtonian cooling which accounts in a crude yet effective way for the emission and the absorption of long and short wave radiation and for the heat convergence associated to convective processes (following Held and Suarez (1994)). This process is described by the equations
(6)  
(7) 
where is the temperature restoration field that depends on the fixed meridional poletoequator temperature gradient and the poletopole gradient . The latter gradient is zero in our experiments, so that we have equatorial symmetry in our boundary conditions and each solution we find is accompanied by a mirrored solution at the equator. For completeness, we also add the full equations of the restoration field, and we refer the reader to Fraedrich et al. (1998) for a more detailed account:
(8)  
(9)  
(10) 
Here, is the height of the tropopause, whereas is the global constant height of the tropopause. is a technical smoothing parameter. Finally, the hyperdiffusion in Eq. 6 is defined as and parametrizes small scale interactions.
PUMA uses spherical harmonics and grid point fields of the prognostic variables. Utilising the Fourier transform along the zonal direction and a Legendre transformation, PUMA computes the linear terms in spectral space and the nonlinear terms in grid point space. The timestepping scheme is a combination of a leapfrog scheme with RobertAsselin filter.
The PUMA User’s Guide includes more details and a complete description of the exact implementation and form of the various forcings (Fraedrich et al., 1998).
Symbol  Description 

temperature  
reference temperature  
temperature deviation from  
relative vorticity  
divergence  
surface pressure pressure  
geopotential  
time  
longitude, latitude  
sigma vertical coordinate  
vertical velocity in system  
vertical velocity in system  
zonal, meridional component of horizontal velocity  
horizontal velocity with components ,  
Coriolis parameter  
diabatic heating rate  
specific heat of dry air at constant pressure  
adiabatic coefficient 
1.2 Maooam
Although the atmospheric dynamics of both models are largely governed by the same processes, MAOOAM differs in many respects from the standalone PUMA model. Most importantly, the atmosphere of MAOOAM features both a mechanical and a thermal coupling with a shallowwater ocean layer, which is absent in PUMA. Furthermore, MAOOAM is a midlatitude model which uses the quasigeostrophic approximation (Charney and Straus, 1980) on a plane (Vallis, 2006), whereas PUMA is a global primitiveequation model, in which the filtering is applied at a much smaller scale. The representation of the dynamical fields differs accordingly, with MAOOAM adopting a Fourier basis, using products of sine and cosine functions that respect the boundary conditions of a zonally periodic atmosphere over a rectangular ocean basin (De Cruz et al., 2016).
The dynamics of MAOOAM’s twolayer atmosphere is described by the quasigeostrophic vorticity equations, expressed in terms of the streamfunction fields at 250 hPa and at 750 hPa as in Charney and Straus (1980),
(11)  
(12) 
in which the vertical velocity can be eliminated by applying the hydrostatic relation and the ideal gas law, as detailed in (De Cruz et al., 2016).
Following Pierini (2011), the equation of motion for the ocean layer is described by
(13) 
The prognostic equations for the atmospheric and oceanic temperature fields, using an energy balance scheme as in Barsugli and Battisti (1998), are
(14)  
(15) 
The quartic terms in these equations are linearised by decomposing the temperature fields around a spatially and temporally constant equilibrium temperature, and , and solving the quartic equation for the equilibrium temperature (Vannitsem et al., 2015).
The thermal wind relation allows one to link the atmospheric temperature anomaly to the baroclinic streamfunction , more specifically . Hence, the remaining independent dynamical fields are the barotropic atmospheric streamfunction field , defined as , the oceanic streamfunction field , and the temperature anomalies and of the atmosphere and the ocean. The other parameters and variables that feature in the MAOOAM model equations are explained in Table 2.
The model equations are nondimensionalized, and the dynamical fields are expanded in a configurable set of Fourier modes. The MAOOAM code computes the coefficients for the resulting set of ordinary differential equations (ODEs) as algebraic formulae of the wavenumbers. These ODEs are then integrated using a fourthorder RungeKutta integration scheme. We refer the reader to De Cruz et al. (2016) for more details on the expansion of the dynamical fields in terms of Fourier modes, the computation of the coefficients, and the tensorial implementation of the ODEs.
In what follows, we will use a shorthand notation that uses the maximum wavenumbers and to specify the model resolution. If the resolution of the ocean and the atmosphere are the same, the model resolution is referred to as x; otherwise, it is denoted as atm. x, oc. x.
Variable (units)  Description 

, (m s)  streamfunction of the ocean, atmosphere 
(Pa s)  vertical velocity in pressure coordinates 
, (K)  temperature of the ocean, atmosphere; 
, (K)  temperature anomaly of the ocean, atmosphere 
Parameter (units)  Description 
meridional to zonal aspect ratio  
(km)  meridional extent of the domain 
(s)  Coriolis parameter at 45 latitude 
(W m K)  heat transfer coefficient at the oceanatmosphere interface 
(s)  friction coefficient at the bottom of the ocean layer 
, (W m)  insolation coefficient of the ocean, atmosphere 
(s)  friction coefficient at oceanatmosphere interface 
(s)  friction coefficient between the atmospheric layers 
(m)  depth of the ocean layer 
(s)  mechanical oceanatmosphere coupling coefficient 
(J kg K)  gas constant of dry air 
(km)  reduced Rossby deformation radius of the ocean 
(kg m)  density of the ocean 
(W m K)  StefanBoltzmann constant 
(m s Pa)  static stability of the atmosphere 
(m s)  Rossby parameter 
, (J m K)  Specific heat capacity of the ocean layer, atmosphere 
, (K)  constant solution for the temperature of the ocean, atmosphere 
greybody atmosphere emissivity 
2 Methodology
2.1 Computation of the Lyapunov exponents
Let us write the evolution laws of the autonomous system presented in Section 1 in a synthetic form
(16) 
where is a vector containing the entire set of relevant variables = such as temperature or wind velocity, projected on a relevant set of modes as described in Section 1. The function is a nonlinear function of the variables and represents a set of parameters.
Let us consider a small perturbation along the trajectory, , generated by model (16), denoted . Provided that this perturbation is sufficiently small (ideally infinitely small), its dynamics can be described by the linearised equation,
(17) 
and a formal solution can be written as
(18) 
where the matrix is referred to as the resolvent matrix. This matrix is responsible for the amplification or contraction of the errors during the time period . In order to get information independent of the initial or final time, a limit should be taken. Oseledets (1968, 2008) demonstrates that provided that the system is ergodic, the following limit exists for almost all initial conditions ,
(19) 
The Lyapunov exponents are then defined as the natural logarithm of the eigenvalues of . These are usually represented in decreasing order and the full set of exponents is called the Lyapunov spectrum. Other definitions are available but will not be discussed here since we do not use them in this study. These can be found in (Eckmann and Ruelle, 1985; Legras and Vautard, 1995), and in a recent work in the context of the coupled oceanatmosphere (Vannitsem and Lucarini, 2016).
The computation of the backward Lyapunov exponents follows the standard algorithm of (Shimada and Nagashima, 1979; Benettin et al., 1980) based on the GramSchmidt orthogonalisation.

An ensemble of perturbation vectors is randomly initialised.

Every time step, the model propagator is computed from the tangent linear model. This is the matrix that quantifies the transition from one model state into that one time step later.

The model is integrated forward in time, and the propagators are accumulated (multiplied) into a matrix .

Every time steps, is evolved with , and GramSchmidt orthogonalised (using a decomposition). The local Lyapunov spectrum is computed from the diagonal of .

Mean and variance of the local Lyapunov exponents are calculated.
The full Lyapunov spectrum of a model allows us to compute some additional interesting properties of its attractor. One of these is the KaplanYorke or Lyapunov dimension , which provides (a lower bound on) the fractal dimension of the attractor, and is defined as (Frederickson et al., 1983):
(20) 
where is the highest index for which the sum of the largest Lyapunov exponents is still strictly positive.
The second important property of the attractor is the KolmogorovSinai or metric entropy , a quantity that describes the rate of growth of the Shannon entropy (Eckmann and Ruelle, 1985; Boffetta et al., 2002), which characterises the quantity of information necessary to locate the solution on its attractor. Its upper bound is the sum of the positive Lyapunov exponents:
(21) 
with the equality proven for a very particular class of systems, known as Axiom A systems. Here the KS entropy will be assumed to be close to the sum of positive exponents, and hence this sum will be referred to as the KS entropy.
2.2 Large deviation laws
Since the Lyapunov exponents are obtained by considering limiting conditions where the initial perturbations are very small and the time span over which the growth or decay rate is very long, they cannot reasonably be used to study predictability outside such conditions. Finitetime Lyapunov exponents (FTLEs) (e.g. Haller, 2000) have been proposed to address such shortcomings, with the caveat that they do not enjoy the extremely beneficial mathematical properties (especially, normindependence) that characterise the Lyapunov exponents.
In this paper, we focus on the FTLEs and their relation to the asymptotic mean LEs. Hence, we are interested in averages over a time of one backward Lyapunov exponent and its statistics. If a dynamical system is an Axiom A system or – invoking the chaotic hypothesis – one of a certain type of non Axiom A systems, these fluctuations for a finite, but large may be described (based on (Schalge et al., 2013; Pazó et al., 2013; Laffargue et al., 2013)) by a large deviation law (Kifer, 1990; Touchette, 2009). For the finitetime backward LEs and for a large , we will verify the following relation for the distribution of the averages:
(22) 
is the rate function, which is independent of . The rate function can be computed directly from this relation as .
If represents a time series, we have to take the autocorrelation into account. The FTLEs for the models under consideration have a nonzero autocorrelation. To account for this, the time series are decomposed into blocks that are decorrelated. For each LE, we find the smallest block size, called the decorrelation time . The time is chosen to be the time lag when the autocorrelation drops below .
2.3 Experimental design: PUMA
We choose a simple setup of PUMA. In this spirit, we also switch off orography. The system is forced via a constant temperature gradient between the equator and the respective poles, as detailed in Section 1.1. We conduct simulations at a horizontal resolution of T42, which amounts to roughly 250 km. In gridpoint space this corresponds to a Gaussian grid with 64 latitudes and 128 longitudes. In the vertical direction we restrict the resolution to 10 sigma levels. The integration scheme uses a time step of one hour.
The objective of our experiments with PUMA is to compute the backward Lyapunov exponents. For this we perform spinup simulations for 30 years from random initial conditions. We then obtain the first 200 Lyapunov exponents using the Benettin algorithm described in Section 2.1. We allow the algorithm to converge for 5 years and finally obtain a time series of 25 years for all LEs. In order to explore two different chaotic regimes with many positive LEs, we perform two experiments with an equatortopole temperature gradient of 50 K and 60 K, respectively (with ).
Note that in order to compute the Lyapunov exponents, it is necessary to construct the tangent linear of PUMA. We generated parts of the code using the program TAF by FastOpt (Giering and Kaminski, 2003).
2.4 Experimental design: MAOOAM
Table 3 lists the values of the physical parameters that are used in the present study. These are selected to lie within the realistic ranges previously derived by Vannitsem and De Cruz (2014), and correspond to the setup used by De Cruz et al. (2016). In addition, we explore different values of the mechanical oceanatmosphere coupling coefficient and the eddy viscosity coefficients and , as well as a range of model resolutions.
Parameter (unit)  Value  Parameter (unit)  Value 

(km)  
(km)  (kg m)  
(s)  (W m K)  
(W m K)  (m s Pa)  
(s)  (m s)  
(W m)  (J m K)  
(W m)  (J m K)  
(s)  (K)  
(s)  (K)  
(m)  
(J kg K) 
All experiments are performed with the same integration parameters. The time step of 0.2 nondimensional time units corresponds to 32.3 minutes in dimensional units. Before calculating the Lyapunov spectrum, a transient run of nondimensional time units is performed, corresponding to 30726 years. Using the tangent linear model of MAOOAM, the backward Lyapunov exponents are then computed using the algorithm described in Section 2.1. In our simulations, the orthogonalisation is performed every time step, i.e. . The Lyapunov spectrum is computed from simulations of 614 years.
The experiments are performed for different resolutions as discussed in Section 1 and for different dissipation schemes as described below.

nodissip

nodissipreducedstress
For this “reducedstress” experiment, the coupling parameter is reduced to .

dissipation
One of the physical processes that was not included in MAOOAM v1.0 (De Cruz et al., 2016) is the kinematic dissipation of energy due to turbulent diffusion, which becomes increasingly important at smaller spatial scales. This process is parametrized as a dissipation term in the prognostic equations for the atmospheric (barotropic) and oceanic streamfunction, that is proportional to the squared Laplacian of the respective streamfunction:
(23) (24) We adopt the values for the parameters and from Van der Avoird et al. (2002), where they are estimated to be
(25) (26) Furthermore, is set to .

dissipationx10
In this experiment, , but the parameters and are set to a higher value:
(27) (28) 
dissipationreducedstress
This experiment has the same parameters as the “dissipation” experiment, except for the coupling parameter which is reduced to .
Note that these idealised experiments do not take into account any dependence of the eddy (or turbulent) viscosity on the truncation scale as usually done in turbulence (e.g. Lesieur, 1990). However, even in the higher resolutions explored so far, we are still far from the scaling regimes for which these dependences may apply. In addition, the values of the eddy viscosity coefficients used are typically valid for a model configuration running at a spatial scale of the order of 100 km (Van der Avoird et al., 2002), which is smaller than the typical truncation used here. For this reason, we have performed a second experiment with a higher eddy viscosity coefficient. The problem of truncation and representation of subgridscale processes is an important open problem in climate modelling that needs careful attention. This matter falls beyond the scope of the present investigation, but forms the subject of a different study in context of the MAOOAM model (Demaeyer and Vannitsem, 2016). Note that in principle the dissipated kinetic energy should become an input to the thermodynamic equations of the system as positive heat source. As discussed in Lucarini and Fraedrich (2009), neglecting this process can have serious dynamical implications on long temporal scales. Additionally, an imperfect representation of this feedback between dynamics and thermodynamics is one of the sources of serious imperfections on the closure of the energy budget in climate models (Lucarini and Ragone, 2011; Lucarini et al., 2014). This issue will be analysed in future investigations.
3 Results
3.1 Puma
Here we present the results for the two different experiments with PUMA, described in Section 2.3, and discuss our findings.
Figure 1 shows the two different Lyapunov spectra obtained in our experiments with PUMA. The averages were computed from a time series of 25 years of daily finitetime LEs. We can estimate the size of the attractor and by that estimate the degrees of freedom inherent to the attractor with the KaplanYorke dimension , as described in Section 2.1. The number of positive exponents and are shown in the legend of Fig. 1. Our findings confirm earlier results using twolayer QG models that suggested an increase of and the number of positive Lyapunov exponents for a higher meridional temperature gradient (Lucarini et al., 2007; Schubert and Lucarini, 2015).
There are two very small exponents since the model setup is zonally symmetric which creates an additional zero mode (see Schubert and Lucarini (2015) for details). Otherwise, there are not many nearzero LEs in PUMA. There is continuity between the time scales that characterise the QG dynamics on the one hand, and faster, smallerscale motions on the other hand. This shows that the usual assumption of a clear timescale separation adopted when applying the filtering to derive the QG equations is, in fact, rather stretched with respect to reality.
Nevertheless, the 50 K spectrum in comparison to the 60 K spectrum has a much smaller slope where the LE are near zero and negative. This may suggest the presence a longer term regime switching behaviour. One potential source for such a regime change is the switching between blocked and nonblocked states of the midlatitudes atmosphere.
We have computed the blocking rate employing the well established TibaldiMolteni Index (Molteni et al., 1988). We indeed find a higher blocking rate for 50 K (%) than for 60 K (%). We would like to explore this connection further in future studies, especially in the direction of studying the location of the CLVs during blocking (Schubert and Lucarini, 2015).
Next, the existence of a large deviation law for the FTLEs is verified, as described in Section 2.2. The decorrelation time is usually between 1 or 3 days. Therefore the rate function is computed for averages of length days.
In Figs. 5, 5, 9, 9, 7, 7, 3 and 3 the results for the rate functions are shown for the fastest growing LE 1 (Figs. 3 and 3), faster decaying LE 150 (Figs. 5 and 5), the positive and negative nearzero LEs 59 and 64 for (Figs. 7 and 7) and nearzero LEs 66 and 71 for (Figs. 9 and 9). Since Eq. 22 is needed to compute the rate function, a long time series is necessary to estimate the distribution reliably. Our intent is to make at least a qualitative assessment of the convergence rate for . The top panels of these figures show the approximation of the respective distributions obtained via kernel smoothing of the distribution of the blockaveraged LEs. The bottom panels show the rate function for different derived from Eq. 22.
We make the following observations. For all LEs the tendency for convergence of the rate function is visible. Also, the rate functions’ shape is approximately parabolic and the estimates of the rate functions converge to the asymptotic with a comparable speed regardless of the value of the corresponding LE.
We interpret these results as another consequence of the nonexisting clearcut timescale separation in a purely atmospheric model like PUMA. This is in opposition to what was speculated in Schubert and Lucarini (2015), where a primitiveequation model was expected to feature a timescale separation visible in the Lyapunov spectrum. Such a timescale separation would have been an a posteriori justification of the filtering by the QG approximation.
We have shown that in a primitiveequation model with a highdimensional phase space of the size of the attractor is small () in comparison. Nevertheless, the unstable subspace can still be regarded as a highdimensional subspace. We also found sound results regarding the existence of a large deviation law independent of the growth rate of the linear perturbations. In hindsight with respect to the findings in MAOOAM this can be explained with the absence of a clear timescale separation.
3.2 Maooam
The Lyapunov analysis is performed on the set of model configurations described in Section 2.4. Let us first evaluate the impact of the resolution on the amplitude of the dominant Lyapunov exponent. The largest Lyapunov exponent , which largely determines the limit of predictability, is plotted as a function of the model resolution for each experiment in Fig. 10. The dominant exponent does not display a clear upward or downward trend versus model resolution and seems to stabilise for higher resolutions. This interesting feature suggests that the lowerorder systems explored here already display a qualitatively correct amplitude for the dominant instability.
Furthermore, as we could expect, the predictability is enhanced for models where the scaledependent dissipation term is present. The decrease in also appears to suggest an enhanced predictability for models which have a larger oceanatmosphere coupling parameter , but this feature is not so clear for higher resolution versions. Vannitsem (2017) studied the dependence of the predictability on this coupling parameter in the loworder 36variable model that lies at the basis of MAOOAM. Two distinct mechanisms were identified to drive the increase in predictability with increasing . To a first approximation, the mechanical coupling of the fast atmosphere to the slow ocean corresponds to an effective friction term which reduces error growth in the atmosphere. Moreover, increasing the oceanatmosphere coupling above a critical value induces a sudden jump in predictability, associated with the development of a slow coupled oceanatmosphere mode (Vannitsem et al., 2015; Vannitsem, 2017).
Figures 11 to 13 show the full sets of Lyapunov exponents, or Lyapunov spectra, for the different experiments. These figures reveal the presence of three ranges in the spectrum of Lyapunov exponents: the positive, negative nearzero and large amplitude negative Lyapunov exponents, associated to the unstable, central and stable manifolds, respectively, in qualitative agreement with what was found in Vannitsem and Lucarini (2016). We expect that the stable and unstable manifolds mainly characterise the dissipative and unstable motions of the atmosphere, while the central manifold also projects considerably on the variables of the ocean.
The highly populated central manifold of MAOOAM is in stark contrast with the few nearzero LEs in PUMA. Being a purely atmospheric model, PUMA’s Lyapunov spectrum does not exhibit the large timescale separation present in MAOOAM. Indeed, the spectrum of PUMA bears more resemblance to that of the QG twolayer model of Schubert (2015).
Upon increasing the number of modes in the ocean and the atmosphere, the number of positive Lyapunov exponents (indicated with a vertical arrow) consistently increases, but not as much as the number of strongly negative exponents. This suggests that most of the additional spatial scales that are resolved by the higherresolution models are highly dissipative, hence increasing the number of strongly negative Lyapunov exponents. The additional positive and nearzero exponents that are introduced at these scales nevertheless indicate that the added resolution still resolves some scales that are important for the description of the dynamics. This is in agreement with the conclusion in De Cruz et al. (2016), where it was shown that in order to describe the ocean dynamics, one needs to be able to resolve the Rhines scale , requiring oceanic wavenumbers as high as of 40–50.
Figure 15 plots the KaplanYorke dimension as a function of the model resolution. This shows that is the highest for the models which do not include the scaledependent dissipation process (nodissip). A reduction in the oceanatmosphere coupling appears to slightly increase for most model resolutions, both in the case with and without scaledependent dissipation. The tenfold increase in the dissipation parameters and (dissipationx10) results in the lowest values for , as can be expected from a more dissipative but still chaotic system.
As the number of dimensions increases quadratically and not linearly for the consecutive model resolutions, it is instructive to rescale by the number of dimensions , as shown in Fig. 16. This shows that while increases with resolution, the attractor dimension’s fraction of the full phase space dimension decreases (even if slowly) with increasing resolution from the atm. 6x6, oc. 6x6 models onward, for all experiments, suggesting that one is adding in higher proportion highly stable modes that do not necessarily play an important role in the dynamics. In other words, we are not in the regime where the system is extensive, as, in fact, the geometry of the domain is fixed and we are capturing a larger and larger (yet insufficient) fraction of the active dynamical processes as the resolution is increased. Had we reached the optimal resolution, Fig. 15 would be flat, and Fig. 16 would approach zero.
Figure 17 shows the KolmogorovSinai entropy versus model resolution, for the different experiments. The trends for the “nodissip” and “nodissipreducedstress” experiments appear to suggest that would increase unboundedly for increasing model resolution if a parametrization for the scaledependent dissipation is absent. The experiments which take this process into account, paint a more realistic picture, with levelling off at the highest model resolutions.
An additional experiment is performed by increasing the resolution of the ocean and of the atmosphere separately starting from a specific symmetric configuration “6x6”. Figure 18 displays the Lyapunov spectra for the model configuration “dissipation”. Two important features stand out: (i) when the resolution of the atmosphere is increased, the majority of the new exponents populate the stable manifold; (ii) on the contrary when the resolution of the ocean is increased, the number of slightly positive and slightly negative exponents increases considerably. This also suggests that the increase of Lyapunov dimension and the number of positive exponents after a resolution of atm. 6x6  oc. 6x6 should be attributed to the presence of the ocean. In this sense the ocean plays an active role in the development of the coupled dynamics. This result deserves an extensive investigation by looking at the properties of the CLVs.
As a final analysis, we have given a preliminary look at whether large deviation laws can be established for the longterm statistics of the FTLEs. In what follows, we consider the “9x9” simulations. Similarly to what was found in a previous analysis performed on a severely truncated version of MAOOAM (Vannitsem and Lucarini, 2016), we find that the time series of the FTLEs corresponding to the strongly damped mode are weakly correlated and one can construct the rate functions defining the large deviations laws; compare figures 19ac) for the 351 LE. Additionally, the lagged time correlation of the nearzero LEs are very strong and it makes no sense to look for large deviation laws in this case.
In contrast to what was presented in Vannitsem and Lucarini (2016), establishing large deviation laws for the FTLEs associated with positive LEs is not trivial, even when one considers the first FTLE. Laggedtime correlations are such that the available time series are not sufficiently long to reach the asymptotic limit, except for the nodissip simulation scenario, which allows for the presence of a larger value of the 1 LE and faster decay of correlations; compare Figs. 20a)c). This suggests that when many unstable modes are present, disentangling their longterm properties requires very long integrations, possibly as a result of geometrical quasidegeneracies among such modes. This is an issue that should be further explored, given its practical and theoretical relevance. One can conjecture that the damped modes do not feature such properties as their dynamics is mostly driven by linear dissipative processes. Therefore, we propose that an accurate analysis of the tangent space with the formalism of CLVs is required to advance our understanding of predictability at medium and long time scales.
In brief, these results indicate that the dominant instabilities of the coupled oceanatmosphere system are well captured by MAOOAM, even at low resolutions. However, the increase of the Lyapunov dimension with the resolution implies that the relevant dynamics of the system are not yet fully resolved, in agreement with De Cruz et al. (2016). The main role of the ocean in this matter is confirmed by varying the ocean and atmosphere resolutions independently. Conversely, increasing the resolution of the atmosphere only adds highly dissipative modes. Finally, in contrast with what was found for a loworder version of MAOOAM (Vannitsem and Lucarini, 2016), large deviation laws cannot be established for the nearzero and positive FTLEs in the “9x9” configuration.
[Toward a new programme] The chaotic nature of the atmosphere and of the climate system has been investigated in the present work in the context of a primitiveequation atmospheric model and a coupled oceanatmosphere model. Both systems suggest that highdimensional dynamical processes are at play with very interesting distinct specificities.
The Lyapunov spectra of the two models considered here have rather different qualitative features, as a result of their structural differences, which have profound impacts on the type of possible instability mechanisms. Following Gallavotti and Lucarini (2014), one expects that if a clear timescale separation between distinct dynamical regimes is present, one should find that the Lyapunov exponents can be divided into separate groups, corresponding to distinct clusters in their values. This is the analogue in full nonlinear terms of what is envisioned by the usual scale analysis of GFD equations.
MAOOAM is a coupled quasigeostrophic atmosphereocean model, which, by definition, features a large timescale separation between ocean and atmosphere, and lacks a satisfactory representation of mesoscale and submesoscale processes. PUMA is an atmosphericonly primitiveequation model, which can represent the faster, smallerscale instabilities associated with processes occurring well below the Rossby deformation radius. On the other side, the lack of an active ocean component removes the presence of very slow scales and does not allow for a builtin scale separation in the dynamics.
We summarise here some findings:

In PUMA the spectrum of Lyapunov exponents changes in accordance with the paradigm that stronger baroclinic forcing leads to a more unstable atmosphere, as already observed in Schubert and Lucarini (2015) for a quasigeostrophic model. The model does not feature any separation of scale, as the Lyapunov spectrum is quite smooth. As a result, one cannot clearly distinguish the modes corresponding to baroclinic instability, KelvinHelmholtz instability, etc. Despite this, we find that for the lower meridional temperature gradient () the spectrum features more negative closer to zero exponents. Interestingly, this might be related to the presence of blocking, but specific studies with more robust dynamics between blocking and nonblocking situations are necessary to clarify that. Additionally, one finds that all the FTLEs accurately obey large deviation laws defining the predictability properties at long time scales, including the nearzero exponents. The model can be categorised as being nonuniformly hyperbolic with a trivial central manifold (the direction of the flow).

For MAOOAM the Lyapunov spectrum is shaped considerably by the presence of the ocean, with a large portion of exponents close to zero. The subspace associated with these exponents corresponds to the central manifold as in the theory of partially hyperbolic systems, and presents features analogous to what was observed in Vannitsem and Lucarini (2016). Furthermore, raising the ocean resolution in MAOOAM clearly increases the number of both positive and negative nearzero Lyapunov exponents, which implies a considerable increase of the Lyapunov dimension of the attractor. Reducing the intensity of the dissipative processes leads, as expected, to an increase in the instability of the model.
One can also conjecture that the set of physical modes, as defined by Yang and Radons (2013), are not yet fully populated since one would expect that the isolated modes are strongly dissipative. This might imply that the resolution necessary to correctly describe the dynamics of the system is much higher in the ocean. This aspect is well known in ocean dynamics since the unstable baroclinic modes, that play an important role in the ocean variability, can only be resolved with scales smaller than 50 km. Yet the question of defining an appropriate resolution (or more exactly an appropriate set of dynamical modes) for which the dynamics is well captured is still open and the analysis of the CLVs in the spirit of (Yang and Radons, 2013; Vannitsem and Lucarini, 2016) can help answer this very important question.
The analysis of the FTLEs of MAOOAM reveals some interesting insight into the dynamics. The FTLEs associated to the strongly dissipative modes obey large deviation laws, while those corresponding to the nearzero LEs do not. This behaviour is expected, and in agreement with what was found in Vannitsem and Lucarini (2016). Surprisingly, however, it is hard to find convergence for the FTLEs associated to the positive LEs. This may point to the presence of nontrivial ocean influence on the (mostly) atmospheric instabilities.
In the programme we want to develop starting from this investigation, we will employ CLVs in highdimensional models to tackle various open problems. CLVs allow to associate growth and decay rates to timedependent physical modes, and provide a geographical portrait of where instability or damping develops.
First, what is the minimal but sufficient resolution? This is a crucial question, in particular in view of the current computer power needed to perform longterm numerical integrations. A possible way to quantify where this threshold might be, is by means of the different modes identified by Yang and Radons (2013) using CLVs. The CLVs provide information on the optimal splitting of physical modes that effectively describe the dynamics of the system and the highly damped modes. The latter can be considered as noisy, purely dissipative terms whose resolution is not necessarily relevant, and are also called isolated modes. Yang and Radons (2013) interpreted them as the result of having a larger number of degrees of freedom in the model than required to resolve all meaningful physical processes. The central feature that allows for the splitting is the angle between the CLVs. If two CLVs display angles around 90 degrees and bounded away from zero degrees, these directions in phase space can be naturally split. Being able to describe the physical modes is deemed essential to satisfactorily reproduce the socalled inertial manifold. The inertial manifold contains the effective finitedimensional dynamics of the system, which, we remind, is originally infinitedimensional if we represent a continuum system described by (S)PDEs. In particular, it would be interesting to determine how the threshold for resolving the inertial manifold varies in a purely atmospheric model compared to a coupled model atmosphereocean model. Additionally, the analysis of geometry of the tangent space can clarify to what extent a system can be treated  effectively, not rigorously  as hyperbolic vs partially hyperbolic. As explained in Vannitsem and Lucarini (2016), this has profound implications for the predictability.
Second, we want to understand multiscale instabilities better and find out what are the driving processes behind their growth. Here, the covariance of the CLVs with the tangent linear equation is the key for understanding instabilities and their properties far away from an equilibrium. Traditionally, even in a chaotic setting such an analysis relied on classic normal mode instability of fixed stationary states (e.g. Charney, 1947; Eady, 1949; Pedlosky, 1964) to explain phenomena like the baroclinic and barotropic instability. This approach has been very beneficial but has many known shortcomings for the understanding of highly nonlinear phenomena such as wavewave interactions (Speranza and Malguzzi, 1988) or regime switching like blocking (Pelly and Hoskins, 2003). Additionally, CLVs will allow us to better understand coupled oceanatmospheric modes. We wish to develop our future programme in line with Schubert and Lucarini (2015, 2016) who demonstrated that CLVs give a picture of what types of instabilities exist in an atmospheric quasigeostrophic (QG) twolayer model and of the energetics behind them. For example, the fastest modes can be almost exclusively barotropically unstable even though traditional normal mode analysis suggests the most unstable modes are driven by baroclinic energy conversion. Given these findings, we expect an even more diverse mixture of different types of instabilities in multiscale systems such as PUMA or MAOOAM. This approach is a promising alternative to restricting the analysis to either studying idealised life cycles of instabilities (Plougonven and Zhang, 2014) or studying yet again normal modes (Molemaker et al., 2005).
The PUMA model is a part of PLASIM, for which the source code can be downloaded at https://www.mi.unihamburg.de/en/arbeitsgruppen/theoretischemeteorologie/modelle/sources/plasim.tgz. The source code for MAOOAM v1.2 is available at http://github.com/Climdyn/MAOOAM. This version is also archived at http://dx.doi.org/10.5281/zenodo.231162. The source code to compute the Lyapunov exponents is available upon request to the corresponding author.
The Lyapunov spectra of the different PUMA and MAOOAM model configurations, that were computed using the Benettin algorithm, are available as supplementary material.
V. Lucarini and S. Schubert performed the analysis of PUMA. S. Schubert wrote the code to compute the LEs for PUMA. L. De Cruz and S. Vannitsem performed the analysis of MAOOAM. J. Demaeyer and S. Schubert wrote the code to compute the LEs for MAOOAM. L. De Cruz, S. Vannitsem and J. Demaeyer wrote MAOOAM. V. Lucarini analysed the FTLEs in MAOOAM. All authors contributed to the writing of the manuscript.
S. Vannitsem and V. Lucarini are members of the editorial board of the journal. The other authors declare that they have no conflict of interest.
Acknowledgements.
LDC was supported by the BELSPO under contract BR/165/A2/Mass2Ant. JD was supported by BELSPO under contract BR/121/A2/STOCHCLIM. VL and SS were supported by the DFG through the contract SFB/Transregio TRR181.References
 Barsugli and Battisti (1998) Barsugli, J. J. and Battisti, D. S.: The Basic Effects of AtmosphereOcean Thermal Coupling on Midlatitude Variability*, Journal of the Atmospheric Sciences, 55, 477–493, https://doi.org/10.1175/15200469(1998)055<0477:TBEOAO>2.0.CO;2, 1998.
 Benettin et al. (1980) Benettin, G., Galgani, L., Giorgilli, A., and Strelcyn, J.M.: Lyapunov Characteristic Exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. Part 1: Theory, Meccanica, 15, 9–20, https://doi.org/10.1007/BF02128236, 1980.
 Boffetta et al. (2002) Boffetta, G., Cencini, M., Falcioni, M., and Vulpiani, A.: Predictability: a way to characterize complexity, Physics Reports, 356, 367–474, https://doi.org/10.1016/S03701573(01)000254, URL https://arxiv.org/pdf/nlin/0101029.pdf, 2002.
 Boschi et al. (2013) Boschi, R., Lucarini, V., and Pascale, S.: Bistability of the climate around the habitable zone: A thermodynamic investigation, Icarus, 226, 1724–1742, https://doi.org/10.1016/j.icarus.2013.03.017, 2013.
 Carrassi et al. (2008) Carrassi, A., Trevisan, A., Descamps, L., Talagrand, O., and Uboldi, F.: Controlling instabilities along a 3DVar analysis cycle by assimilating in the unstable subspace: a comparison with the EnKF, Nonlinear Processes in Geophysics, 15, 503–521, https://doi.org/10.5194/npg155032008, 2008.
 Cencini and Ginelli (2013) Cencini, M. and Ginelli, F.: Lyapunov analysis: from dynamical systems theory to applications, Journal of Physics A: Mathematical and Theoretical, 46, 250 301, https://doi.org/10.1088/17518113/46/25/250301, 2013.
 Cencini et al. (2010) Cencini, M., Cecconi, F., and Vulpiani, A.: Chaos: From Simple Models to Complex Systems, vol. 17 of Series on Advances in Statistical Mechanics, World Scientific, https://doi.org/10.1142/7351, 2010.
 Charney (1947) Charney, J. G.: The Dynamics Of Long Waves In A Baroclinic Current, Journal of Meteorology, 4, 136–162, https://doi.org/10.1175/15200469(1947)004<0136:TDOLWI>2.0.CO;2, 1947.
 Charney and Straus (1980) Charney, J. G. and Straus, D. M.: FormDrag Instability, Multiple Equilibria and Propagating Planetary Waves in Baroclinic, Orographically Forced, Planetary Wave Systems, Journal of the Atmospheric Sciences, 37, 1157–1176, https://doi.org/10.1175/15200469(1980)037<1157:FDIMEA>2.0.CO;2, 1980.
 De Cruz et al. (2016) De Cruz, L., Demaeyer, J., and Vannitsem, S.: The Modular ArbitraryOrder OceanAtmosphere Model: MAOOAM v1.0, Geoscientific Model Development, 9, 2793–2808, https://doi.org/10.5194/gmd927932016, 2016.
 Demaeyer and Vannitsem (2016) Demaeyer, J. and Vannitsem, S.: Stochastic parametrization of subgridscale processes in coupled ocean–atmosphere systems: benefits and limitations of response theory, Quarterly Journal of the Royal Meteorological Society, 143, 881–896, https://doi.org/10.1002/qj.2973, 2016.
 Eady (1949) Eady, E. T.: Long Waves and Cyclone Waves, Tellus, 1, 33–52, https://doi.org/10.1111/j.21533490.1949.tb01265.x, 1949.
 Eckmann and Ruelle (1985) Eckmann, J. and Ruelle, D.: Ergodic theory of chaos and strange attractors, Reviews of Modern Physics, 57, 617–656, https://doi.org/10.1103/RevModPhys.57.617, 1985.
 Fraedrich and Kirk (2005) Fraedrich, K. and Kirk, E.: The portable university model of the atmosphere (PUMA): Storm track dynamics and lowfrequency variability, Meteorologische Zeitschrift, 14, 735–745, https://doi.org/10.1127/09412948/2005/0074, 2005.
 Fraedrich et al. (1998) Fraedrich, K., Kirk, E., and Lunkeit, F.: Portable University Model of the Atmosphere, Tech. rep., DKRZ, Hamburg, URL https://www.dkrz.de/mms//pdf/reports/ReportNo.16.pdf, 1998.
 Frederickson et al. (1983) Frederickson, P., Kaplan, J. L., Yorke, E. D., and Yorke, J. A.: The Liapunov dimension of strange attractors, Journal of Differential Equations, 49, 185–207, https://doi.org/10.1016/00220396(83)900116, 1983.
 Gallavotti and Cohen (1995) Gallavotti, G. and Cohen, E. G. D.: Dynamical Ensembles in Nonequilibrium Statistical Mechanics, Physical Review Letters, 74, 2694–2697, https://doi.org/10.1103/PhysRevLett.74.2694, 1995.
 Gallavotti and Lucarini (2014) Gallavotti, G. and Lucarini, V.: Equivalence of Nonequilibrium Ensembles and Representation of Friction in Turbulent Flows: The Lorenz 96 Model, Journal of Statistical Physics, 156, 1027–1065, https://doi.org/10.1007/s1095501410516, URL http://arxiv.org/abs/1404.6638, 2014.
 Gallez and Babloyantz (1991) Gallez, D. and Babloyantz, A.: Lyapunov exponents for nonuniform attractors, Physics Letters, 161, 247–254, https://doi.org/10.1016/03759601(91)90012W, 1991.
 Giering and Kaminski (2003) Giering, R. and Kaminski, T.: Applying TAF to generate efficient derivative code of Fortran 7795 programs, PAMM, 2, 54–57, https://doi.org/10.1002/pamm.200310014, 2003.
 Ginelli et al. (2007) Ginelli, F., Poggi, P., Turchi, A., Chaté, H., Livi, R., and Politi, A.: Characterizing dynamics with covariant Lyapunov vectors, Physical Review Letters, 99, 130 601, https://doi.org/10.1103/PhysRevLett.99.130601, URL http://arxiv.org/abs/0706.0510, 2007.
 Haller (2000) Haller, G.: Finding finitetime invariant manifolds in twodimensional velocity fields, Chaos: An Interdisciplinary Journal of Nonlinear Science, 10, 99–108, 2000.
 Held and Suarez (1994) Held, I. M. and Suarez, M. J.: A Proposal for the Intercomparison of the Dynamical Cores of Atmospheric General Circulation Models, Bulletin of the American Meteorological Society, 75, 1825–1830, https://doi.org/10.1175/15200477(1994)075<1825:APFTIO>2.0.CO;2, 1994.
 Holden et al. (2013) Holden, P. B., Edwards, N. R., Garthwaite, P. H., Fraedrich, K., Lunkeit, F., Kirk, E., Labriet, M., Kanudia, A., and Babonneau, F.: PLASIMENTSem v1.0: a spatiotemporal emulator of future climate change for impacts assessment, Geoscientific Model Development, 7, 433–451, https://doi.org/10.5194/gmd74332014, 2013.
 Holton (2004) Holton, J. R.: An Introduction To Dynamic Meteorology, Academic Press, fourth edn., 2004.
 Kalnay (2003) Kalnay, E.: Atmospheric Modeling, Data Assimilation, and Predictability, Cambridge University Press, 2003.
 Kifer (1990) Kifer, Y.: Large deviations in dynamical systems and stochastic processes, Transactions of the American Mathematical Society, 321, 505–524, https://doi.org/10.1090/S00029947199010257567, 1990.
 Klein (2010) Klein, R.: ScaleDependent Models for Atmospheric Flows, Annual Review of Fluid Mechanics, 42, 249–274, https://doi.org/10.1146/annurevfluid121108145537, 2010.
 Laffargue et al. (2013) Laffargue, T., Lam, K.D. N. T., Kurchan, J., and Tailleur, J.: Large deviations of Lyapunov exponents, Journal of Physics A: Mathematical and Theoretical, 46, 254 002, https://doi.org/10.1088/17518113/46/25/254002, 2013.
 Legras and Ghil (1985) Legras, B. and Ghil, M.: Persistent Anomalies, Blocking and Variations in Atmospheric Predictability, Journal of the Atmospheric Sciences, 42, 433–471, https://doi.org/10.1175/15200469(1985)042<0433:PABAVI>2.0.CO;2, 1985.
 Legras and Vautard (1995) Legras, B. and Vautard, R.: A Guide to Liapunov vectors, in: Proceedings 1995 ECMWF Seminar on Predictability, vol. 1, pp. 143–156, URL http://gershwin.ens.fr/legras/publis/liapunov.pdf, 1995.
 Lesieur (1990) Lesieur, M.: Introduction to Turbulence in Fluid Mechanics, in: Turbulence in Fluids: Stochastic and Numerical Modelling, pp. 1–18, Springer, Dordrecht, https://doi.org/10.1007/9789400905337_1, 1990.
 Lorenz (1963) Lorenz, E. N.: Deterministic Nonperiodic Flow, Journal of the Atmospheric Sciences, 20, 130–141, https://doi.org/10.1175/15200469(1963)020<0130:DNF>2.0.CO;2, 1963.
 Lucarini and Fraedrich (2009) Lucarini, V. and Fraedrich, K.: Symmetry breaking, mixing, instability, and lowfrequency variability in a minimal Lorenzlike system, Physical Review E, 80, 026 313, https://doi.org/10.1103/PhysRevE.80.026313, 2009.
 Lucarini and Ragone (2011) Lucarini, V. and Ragone, F.: ENERGETICS OF CLIMATE MODELS: NET ENERGY BALANCE AND MERIDIONAL ENTHALPY TRANSPORT, Reviews of Geophysics, 49, RG1001, https://doi.org/10.1029/2009RG000323, 2011.
 Lucarini et al. (2007) Lucarini, V., Speranza, A., and Vitolo, R.: Parametric smoothness and selfscaling of the statistical properties of a minimal climate model: What beyond the mean field theories?, Physica D: Nonlinear Phenomena, 234, 105–123, https://doi.org/10.1016/j.physd.2007.07.006, 2007.
 Lucarini et al. (2014) Lucarini, V., Blender, R., Herbert, C., Ragone, F., Pascale, S., and Wouters, J.: Mathematical and physical ideas for climate science, Reviews of Geophysics, 52, 809–859, https://doi.org/10.1002/2013RG000446, URL http://arxiv.org/abs/1311.1190, 2014.
 Lucarini et al. (2017) Lucarini, V., Ragone, F., and Lunkeit, F.: Predicting climate change using response theory: global averages and spatial patterns, Journal of Statistical Physics, 166, 1036–1064, https://doi.org/10.1007/s109550161506z, 2017.
 Manneville (1985) Manneville, P.: Liapounov exponents for the KuramotoSivashinsky model, in: Macroscopic Modelling of Turbulent Flows, vol. 230, pp. 319–326, Springer Berlin Heidelberg, Berlin, Heidelberg, https://doi.org/10.1007/3540156445_26, 1985.
 Manneville (1995) Manneville, P.: Dissipative structures and weak turbulence, in: Chaos — The Interplay Between Stochastic and Deterministic Behaviour, vol. 457, pp. 257–272, Springer Berlin Heidelberg, Berlin, Heidelberg, https://doi.org/10.1007/3540601880_59, 1995.
 Molemaker et al. (2005) Molemaker, M. J., McWilliams, J. C., and Yavneh, I.: Baroclinic Instability and Loss of Balance, Journal of Physical Oceanography, 35, 1505–1517, https://doi.org/10.1175/JPO2770.1, 2005.
 Molteni et al. (1988) Molteni, F., Sutera, A., and Tronci, N.: The EOFs of the Geopotential Eddies at 500 mb in Winter and Their Probability Density Distributions, Journal of the Atmospheric Sciences, 45, 3063–3080, https://doi.org/10.1175/15200469(1988)045<3063:TEOTGE>2.0.CO;2, 1988.
 Nese and Dutton (1993) Nese, J. M. and Dutton, J. A.: Quantifying predictability variations in a loworder occanatmosphere model: a dynamical systems approach, Journal of Climate, 6, 185–204, https://doi.org/10.1175/15200442(1993)006<0185:QPVIAL>2.0.CO;2, 1993.
 Nicolis (2007) Nicolis, C.: Dynamics of model error: The role of the boundary conditions, Journal of the Atmospheric Sciences, 64, 204–215, https://doi.org/doi.org/10.1175/JAS3806.1, 2007.
 Nicolis et al. (1992) Nicolis, C., Nicolis, G., and Wang, Q.: SENSITIVITY TO INITIAL CONDITIONS IN SPATIALLY DISTRIBUTED SYSTEMS, International Journal of Bifurcation and Chaos, 02, 263–269, https://doi.org/10.1142/S0218127492000276, 1992.
 Nicolis et al. (2009) Nicolis, C., Perdigao, R. A., and Vannitsem, S.: Dynamics of prediction errors under the combined effect of initial condition and model errors, Journal of the Atmospheric Sciences, 66, 766–778, https://doi.org/10.1175/2008JAS2781.1, 2009.
 Oseledets (2008) Oseledets, V.: Oseledets theorem, Scholarpedia, 3, 1846, https://doi.org/10.4249/scholarpedia.1846, 2008.
 Oseledets (1968) Oseledets, V. I.: A multiplicative ergodic theorem. Characteristic Ljapunov, exponents of dynamical systems, Transactions of the Moscow Mathematical Society, 19, 179–210, URL http://www.mathnet.ru/php/archive.phtml?wshow=paper&jrnid=mmo&paperid=214&option_lang=eng, 1968.
 Ott (2002) Ott, E.: Chaos in Dynamical Systems  2nd Edition, Cambridge University Press, https://doi.org/10.1017/CBO9780511803260, 2002.
 Pazó et al. (2008) Pazó, D., Szendro, I. G., López, J. M., and Rodríguez, M. a.: Structure of characteristic Lyapunov vectors in spatiotemporal chaos, Physical Review E, 78, 1–9, https://doi.org/10.1103/PhysRevE.78.016209, URL http://arxiv.org/abs/0808.0432, 2008.
 Pazó et al. (2010) Pazó, D., Rodríguez, M. A., and López, J. M.: Spatiotemporal evolution of perturbations in ensembles initialized by bred, Lyapunov and singular vectors, Tellus, Series A: Dynamic Meteorology and Oceanography, 62, 10–23, https://doi.org/10.1111/j.16000870.2009.00419.x, 2010.
 Pazó et al. (2013) Pazó, D., López, J. M., and Politi, A.: Universal scaling of Lyapunovexponent fluctuations in spacetime chaos, Physical Review E, 87, 062 909, https://doi.org/10.1103/PhysRevE.87.062909, 2013.
 Pedlosky (1964) Pedlosky, J.: The Stability of Currents in the Atmosphere and the Ocean: Part II, Journal of the Atmospheric Sciences, 21, 201–219, https://doi.org/10.1175/15200469(1964)021{%}3C0342{%}3ATSOCIT{%}3E2.0.CO{%}3B2, 1964.
 Pelly and Hoskins (2003) Pelly, J. L. and Hoskins, B. J.: A New Perspective on Blocking, Journal of the Atmospheric Sciences, 60, 743–755, https://doi.org/10.1175/15200469(2003)060<0743:ANPOB>2.0.CO;2, 2003.
 Pesin (2004) Pesin, Y.: Lectures on Partial Hyperbolicity and Stable Ergodicity, European Mathematical Society, 2004.
 Pierini (2011) Pierini, S.: Lowfrequency variability, coherence resonance, and phase selection in a loworder model of the winddriven ocean circulation, Journal of Physical Oceanography, 41, 1585–1604, https://doi.org/10.1175/JPOD1005018.1, 2011.
 Plougonven and Zhang (2014) Plougonven, R. and Zhang, F.: Internal gravity waves from atmospheric jets and fronts, Reviews of Geophysics, 52, 33–76, https://doi.org/10.1002/2012RG000419, 2014.
 Ragone et al. (2015) Ragone, F., Lucarini, V., and Lunkeit, F.: A new framework for climate sensitivity and prediction: a modelling perspective, Climate Dynamics, pp. 1–13, https://doi.org/10.1007/s0038201526573, 2015.
 Ruelle (1979) Ruelle, D.: Ergodic theory of differentiable dynamical systems, Publications mathématiques de l’IHÉS, 50, 27–58, https://doi.org/10.1007/BF02684768, 1979.
 Ruelle (2009) Ruelle, D.: A review of linear response theory for general differentiable dynamical systems, Nonlinearity, 22, 855–870, https://doi.org/10.1088/09517715/22/4/009, 2009.
 Schalge et al. (2013) Schalge, B., Blender, R., Wouters, J., Fraedrich, K., and Lunkeit, F.: Evidence for a fluctuation theorem in an atmospheric circulation model, Physical Review E, 8760, https://doi.org/10.1103/PhysRevE.87.052113, URL https://mipub.cen.unihamburg.de/fileadmin/files/forschung/theomet/docs/pdf_2013/2013SchalgeflucttheorPRE.pdf, 2013.
 Schneider (2006) Schneider, T.: The General Circulation of the Atmosphere, Annual Review of Earth and Planetary Sciences, 34, 655–688, https://doi.org/10.1146/annurev.earth.34.031405.125144, 2006.
 Schubert (2015) Schubert, S.: Statistical Mechanics of the Fluctuations of a Turbulent QuasiGeostrophic Model of the Atmosphere : Instabilities and Feedbacks, Ph.D. thesis, University Of Hamburg, https://doi.org/urn:nbn:de:gbv:1877942, URL http://ediss.sub.unihamburg.de/volltexte/2016/7794/, 2015.
 Schubert and Lucarini (2015) Schubert, S. and Lucarini, V.: Covariant Lyapunov vectors of a quasigeostrophic baroclinic model: Analysis of instabilities and feedbacks, Quarterly Journal of the Royal Meteorological Society, 141, 3040–3055, https://doi.org/10.1002/qj.2588, 2015.
 Schubert and Lucarini (2016) Schubert, S. and Lucarini, V.: Dynamical analysis of blocking events: spatial and temporal fluctuations of covariant Lyapunov vectors, Quarterly Journal of the Royal Meteorological Society, 142, 2143–2158, https://doi.org/10.1002/qj.2808, 2016.
 Shimada and Nagashima (1979) Shimada, I. and Nagashima, T.: A Numerical Approach to Ergodic Problem of Dissipative Dynamical Systems, Progress of Theoretical Physics, 61, 1605–1616, https://doi.org/10.1143/PTP.61.1605, 1979.
 Speranza and Malguzzi (1988) Speranza, A. and Malguzzi, P.: The Statistical Properties of a Zonal Jet in a Baroclinic Atmosphere: A Semilinear Approach. Part I: Quasigeostrophic, TwoLayer Model Atmosphere, Journal of the Atmospheric Sciences, 45, 3046–3062, https://doi.org/10.1175/15200469(1988)045<3046:TSPOAZ>2.0.CO;2, 1988.
 Sprott (2010) Sprott, J. C.: Elegant Chaos, World Scientific, https://doi.org/10.1142/7183, 2010.
 Thompson (1957) Thompson, P. D.: Uncertainty of initial state as a factor in the predictability of large scale atmospheric flow patterns, Tellus, 9, 275–295, https://doi.org/10.3402/tellusa.v9i3.9111, 1957.
 Touchette (2009) Touchette, H.: The large deviation approach to statistical mechanics, Physics Reports, 478, 1–69, https://doi.org/10.1016/j.physrep.2009.05.002, 2009.
 Trevisan and Pancotti (1998) Trevisan, A. and Pancotti, F.: Periodic Orbits, Lyapunov Vectors, and Singular Vectors in the Lorenz System, Journal of the Atmospheric Sciences, 55, 390–398, https://doi.org/10.1175/15200469(1998)055<0390:POLVAS>2.0.CO;2, 1998.
 Trevisan et al. (2010) Trevisan, A., D’Isidoro, M., and Talagrand, O.: Fourdimensional variational assimilation in the unstable subspace and the optimal subspace dimension, Quarterly Journal of the Royal Meteorological Society, 136, 487–496, 2010.
 Vallis (2006) Vallis, G. K.: Atmospheric and oceanic fluid dynamics : fundamentals and largescale circulation, Cambridge University Press, 2006.
 Van der Avoird et al. (2002) Van der Avoird, E., Dijkstra, H., Nauw, J., and Schuurmans, C.: Nonlinearly induced lowfrequency variability in a midlatitude coupled oceanatmosphere model of intermediate complexity, Climate Dynamics, 19, 303–320, https://doi.org/10.1007/s003820010220x, 2002.
 Vannitsem (2017) Vannitsem, S.: Predictability of largescale atmospheric motions: Lyapunov exponents and error dynamics, Chaos: An Interdisciplinary Journal of Nonlinear Science, 27, 32 101, https://doi.org/10.1063/1.4979042, 2017.
 Vannitsem and De Cruz (2014) Vannitsem, S. and De Cruz, L.: A 24variable loworder coupled ocean–atmosphere model: OAQGWS v2, Geoscientific Model Development, 7, 649–662, https://doi.org/10.5194/gmd76492014, 2014.
 Vannitsem and Lucarini (2016) Vannitsem, S. and Lucarini, V.: Statistical and Dynamical Properties of Covariant Lyapunov Vectors in a Coupled AtmosphereOcean Model  Multiscale Effects, Geometric Degeneracy, and Error Dynamics, Journal of Physics A: Mathematical and Theoretical, 49, 224 001, https://doi.org/10.1088/17518113/49/22/224001, URL http://arxiv.org/abs/1510.00298, 2016.
 Vannitsem and Nicolis (1994) Vannitsem, S. and Nicolis, C.: Predictability experiments on a simplified thermal convection model: The role of spatial scales, Journal of Geophysical Research, 99, 10 377, https://doi.org/10.1029/94JD00248, 1994.
 Vannitsem and Nicolis (1996) Vannitsem, S. and Nicolis, C.: ERROR GROWTH DYNAMICS IN SPATIALLY EXTENDED SYSTEMS, International Journal of Bifurcation and Chaos, 06, 2223–2235, https://doi.org/10.1142/S0218127496001466, 1996.
 Vannitsem and Nicolis (1997) Vannitsem, S. and Nicolis, C.: Lyapunov Vectors and Error Growth Patterns in a T21L3 Quasigeostrophic Model, Journal of the Atmospheric Sciences, 54, 347–361, https://doi.org/10.1175/15200469(1997)054<0347:LVAEGP>2.0.CO;2, 1997.
 Vannitsem et al. (2015) Vannitsem, S., Demaeyer, J., De Cruz, L., and Ghil, M.: Lowfrequency variability and heat transport in a loworder nonlinear coupled oceanatmosphere model, Physica D: Nonlinear Phenomena, 309, 71–85, https://doi.org/10.1016/j.physd.2015.07.006, 2015.
 Wolfe and Samelson (2007) Wolfe, C. L. and Samelson, R. M.: An efficient method for recovering Lyapunov vectors from singular vectors, Tellus A, 59, 355–366, https://doi.org/10.1111/j.16000870.2007.00234.x, 2007.
 Xu and Paul (2016) Xu, M. and Paul, M. R.: Covariant Lyapunov vectors of chaotic RayleighBénard convection, Physical Review E, 93, 062 208, https://doi.org/10.1103/PhysRevE.93.062208, 2016.
 Yamada and Ohkitani (1988) Yamada, M. and Ohkitani, K.: Lyapunov spectrum of a model of twodimensional turbulence, Physical Review Letters, 60, 983–986, https://doi.org/10.1103/PhysRevLett.60.983, 1988.
 Yang and Radons (2013) Yang, H.l. and Radons, G.: Hydrodynamic Lyapunov modes and effective degrees of freedom of extended systems, Journal of Physics A: Mathematical and Theoretical, 46, 254 015, https://doi.org/10.1088/17518113/46/25/254015, 2013.