[1]LesleyDe Cruz [2]SebastianSchubert [1]JonathanDemaeyer
[2,3,4]ValerioLucarini [1]StéphaneVannitsem

1]Royal Meteorological Institute of Belgium, Brussels, Belgium
2]Meteorological Institute, CEN, University of Hamburg, Hamburg, Germany
3]Department of Mathematics and Statistics, University of Reading, Reading, UK
4]Centre for the Mathematics of Planet Earth, University of Reading, Reading,
UK

Lesley De Cruz<EMAIL_ADDRESS>
# Exploring the Lyapunov instability properties of high-dimensional
atmospheric and climate models

###### Abstract

The stability properties of intermediate-order climate models are investigated
by computing their Lyapunov exponents (LEs). The two models considered are
PUMA (Portable University Model of the Atmosphere), a primitive-equation
simple general circulation model, and MAOOAM (Modular Arbitrary-Order Ocean-
Atmosphere Model), a quasi-geostrophic coupled ocean-atmosphere model on a
$\beta$-plane. We wish 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 Kaplan-Yorke dimension of
the attractor increases as well. The convergence rate of the rate function for
the large deviation law of the finite-time Lyapunov exponents (FTLEs) is fast
for all exponents, which can be interpreted as resulting from the absence of a
clear-cut atmospheric time-scale 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 atmosphere-ocean coupling
coefficient or introducing a turbulent diffusion parametrization reduces the
Kaplan-Yorke dimension and Kolmogorov-Sinai entropy. In all considered
configurations, we are not yet in the regime in which one can robustly define
large deviations laws describing the statistics of the FTLEs.

This paper highlights the need to investigate the natural variability of the
atmosphere-ocean 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.

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 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 low-order 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
intermediate-order atmospheric quasi-geostrophic 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 high-dimensional 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 low-order system.

However, 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, 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 surface oceanic circulation, by contrast, is mostly
mechanically driven by atmospheric winds (Lucarini et al., 2014). On even
longer timescales, buoyancy fluxes are an important driver for the deep
ocean’s thermohaline circulation.

This raises the question as to the impact of the coupling to other sub-domains
of the climate system: are the other sub-domains of the climate system
stabilising the atmosphere or not? Vannitsem et al. (2015) partly addressed
this question in the context of coupled low-order ocean-atmosphere 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 full-scale climate system including the different sub-domains is still
open. Recently a new coupled ocean-atmosphere 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
Arbitrary-Order Ocean-Atmosphere 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 co-varies
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 high-dimensional systems. The CLVs have been used to study e.g. spatio-
temporal chaos (Pazó et al., 2008, 2010), Rayleigh-Bénard convection (Xu and
Paul, 2016), and the dynamics of the mid-latitudes atmosphere in the quasi-
geostrophic (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) and allow for associating unstable/stable modes to specific
paths of energy exchange and conversion. 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 and neutral 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
$\pi$ in the ideal case of uniform hyperbolicity. This point of view
complements the investigation of the statistical properties of finite-time
Lyapunov exponents (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
high-dimensional 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 (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 sub-domains with vastly different time scales as done in the case of a
low-order coupled ocean-atmosphere 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 so-called
partially hyperbolic systems (Pesin, 2004). The CLVs corresponding to the
central manifold are geometrically quasi-degenerate, 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 super-exponential behaviour. The question
is therefore what should be the resolution of the coupled atmosphere-ocean
model and what should be the time of observation such that a better separation
emerges between such modes.

###  0.3 Programmatic goals

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
time-dependent perturbations.

In the present manuscript, we explore for the first time the Lyapunov spectra
of the primitive-equation model PUMA, and of intermediate-order configurations
of the coupled ocean-atmosphere 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 Kaplan-Yorke
dimension and Kolmogorov-Sinai 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 Fig. 21 we present our
future programme, which aims to clarify the instabilities of high-resolution
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 state-of-the-art 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), study tipping points (Lucarini
et al., 2010; Boschi et al., 2013) or create large ensembles for climate
change experiments allowing for interesting applications in economic modelling
(Holden et al., 2014) or to assess the feasibility of linear response theory
(Ragone et al., 2016).

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 Navier-Stokes 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 $\eta$, the
divergence $D$, the logarithm of the surface pressure $\ln p_{s}$ and the
temperature $T$. These equations are

$\displaystyle\noindent\begin{split}\frac{\partial\left(\eta+f\right)}{\partial
t}=\frac{1}{1-\mu^{2}}\frac{\partial F_{v}}{\partial\lambda}-\frac{\partial
F_{u}}{\partial\mu}+P_{\eta}\end{split}$ (1)
$\displaystyle\begin{split}\frac{\partial D}{\partial
t}=\frac{1}{1-\mu^{2}}\frac{\partial F_{u}}{\partial\lambda}-\frac{\partial
F_{v}}{\partial\mu}-\left(\frac{U^{2}+V^{2}}{2(1-\mu^{2})}+\Phi+T_{0}\ln
p_{s}\right)+P_{D}\end{split}$ (2) $\displaystyle\begin{split}\frac{\partial
p_{s}}{\partial t}=-\int_{0}^{1}Ad\sigma\end{split}$ (3)
$\displaystyle\begin{split}\frac{\partial T^{\prime}}{\partial
t}=-\frac{1}{1-\mu^{2}}\frac{\partial(UT^{\prime})}{\partial\lambda}-\frac{\partial(VT^{\prime})}{\partial\mu}+DT^{\prime}-\dot{\sigma}\frac{\partial
T}{\partial\sigma}+\kappa\frac{T}{p}\omega+\frac{J}{c_{p}}+P_{T}\end{split}$
(4)

where the vorticity is defined as $\eta=\partial_{x}v-\partial_{y}u$ and the
divergence is defined as $D=\partial_{x}u+\partial_{y}v$. Additionally, one
takes into account the hydrostatic relation

$\frac{\partial\Phi}{\partial\ln\sigma}=-T$ (5)

and $T^{\prime}$ is defined as $T^{\prime}=T-T_{0}$ with $T_{0}=250K$. Some
abbreviations have been used:

$\displaystyle F_{u}$
$\displaystyle=V\left(\eta+f\right)-\dot{\sigma}\frac{\partial
U}{\partial\sigma}-T^{\prime}\frac{\partial\ln p_{s}}{\partial\lambda}$
$\displaystyle F_{v}$
$\displaystyle=-U\left(\eta+f\right)-\dot{\sigma}\frac{\partial
V}{\partial\sigma}-T^{\prime}\left(1-\mu^{2}\right)\frac{\partial\ln
p_{s}}{\partial\mu}$ $\displaystyle A$ $\displaystyle=D+\vec{V}\cdot\nabla\ln
p_{s}$

with $U=u\cos\phi\text{ and }V=v\cos\phi$. 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

$\displaystyle\frac{J}{c_{p}}+P_{T}$
$\displaystyle=\frac{T_{R}\left(\phi,\sigma\right)-T}{\tau_{R}}+H_{T}$ (6)
$\displaystyle T_{R}\left(\phi,\sigma\right)$
$\displaystyle=T^{vert}_{R}(\sigma)+f(\sigma)T^{hor}_{R}(\phi),$ (7)

where $T_{R}$ is the temperature restoration field that depends on the fixed
meridional pole-to-equator temperature gradient $\Delta T_{EP}$ and the pole-
to-pole gradient $\Delta T_{NS}$ . 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:

$\displaystyle T^{hor}_{R}(\phi)$ $\displaystyle=\Delta
T_{NS}\frac{\sin(\phi)}{2}-\Delta T_{EP}\left(\sin^{2}\phi-\frac{1}{3}\right)$
(8) $\displaystyle T^{vert}_{R}(\sigma)$ $\displaystyle=\Delta
T_{trop}+\sqrt{\frac{L}{2}\left(z_{tp}-z(\sigma)\right)^{2}+S^{2}}+\frac{L}{2}\left(z_{tp}-z(\sigma)\right)$
(9) $\displaystyle f(\sigma)$
$\displaystyle=\begin{cases}\sin\left(\frac{\pi}{2}\left(\frac{\sigma-\sigma_{tp}}{1-\sigma_{tp}}\right)\right)&\text{,
if}\ \sigma\geq\sigma_{tp}\\\ 0&\text{, if}\ \sigma<\sigma_{tp}.\end{cases}$
(10)

Here, $\sigma_{tp}$ is the height of the tropopause, whereas $z_{tp}$ is the
global constant height of the tropopause. $S$ is a technical smoothing
parameter. Finally, the hyperdiffusion $H_{T}$ in Eq. 6 is defined as
$H_{T}=\nabla^{8}T$ 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 time-stepping scheme is a
combination of a leap-frog scheme with Robert-Asselin 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).

Table 1: Symbols and variables in the PUMA equations. Symbol | Description
---|---
$T$ | temperature
$T_{0}$ | reference temperature
$T^{\prime}=T-T_{0}$ | temperature deviation from $T_{0}$
$\eta$ | relative vorticity
$D$ | divergence
$p_{s}$ | surface pressure
$\Phi$ | geopotential
$t$ | time
$\lambda,\phi$ | longitude, latitude
$\mu=\sin\phi$ |
$\sigma=\frac{p}{p_{s}}$ | sigma vertical coordinate
$\dot{\sigma}=\frac{d\sigma}{dt}$ | vertical velocity in $\sigma$-system
$\dot{p}=\frac{dp}{dt}$ | vertical velocity in $p$-system
$u,v$ | zonal, meridional component of horizontal velocity
$\vec{V}$ | horizontal velocity with components $U$, $V$
$f$ | Coriolis parameter
$J$ | diabatic heating rate
$c_{p}$ | specific heat of dry air at constant pressure
$\kappa$ | 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 stand-alone PUMA
model. Most importantly, the atmosphere of MAOOAM features both a mechanical
and a thermal coupling with a shallow-water ocean layer, which is absent in
PUMA. Furthermore, MAOOAM is a mid-latitude model which uses the quasi-
geostrophic approximation (Charney and Straus, 1980) on a $\beta$-plane
(Vallis, 2006), whereas PUMA is a global primitive-equation 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 two-layer atmosphere is described by the quasi-
geostrophic vorticity equations, expressed in terms of the streamfunction
fields $\psi^{1}_{\text{a}}$ at 250 hPa and $\psi^{3}_{\text{a}}$ at 750 hPa
as in Charney and Straus (1980),

$\displaystyle\frac{\partial}{\partial
t}\left(\nabla^{2}\psi^{1}_{\text{a}}\right)+J(\psi^{1}_{\text{a}},\nabla^{2}\psi^{1}_{\text{a}})+\beta\frac{\partial\psi^{1}_{\text{a}}}{\partial
x}$
$\displaystyle=-k^{\prime}_{d}\nabla^{2}(\psi^{1}_{\text{a}}-\psi^{3}_{\text{a}})+\frac{f_{0}}{\Delta
p}\omega,$ (11) $\displaystyle\frac{\partial}{\partial
t}\left(\nabla^{2}\psi^{3}_{\text{a}}\right)+J(\psi^{3}_{\text{a}},\nabla^{2}\psi^{3}_{\text{a}})+\beta\frac{\partial\psi^{3}_{\text{a}}}{\partial
x}$
$\displaystyle=+k^{\prime}_{d}\nabla^{2}(\psi^{1}_{\text{a}}-\psi^{3}_{\text{a}})-\frac{f_{0}}{\Delta
p}\omega-k_{d}\nabla^{2}(\psi^{3}_{\text{a}}-\psi_{\text{o}}),$ (12)

in which the vertical velocity $\omega$ 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

$\displaystyle\frac{\partial}{\partial
t}\left(\nabla^{2}\psi_{\text{o}}-\frac{\psi_{\text{o}}}{L_{\text{R}}^{2}}\right)+J(\psi_{\text{o}},\nabla^{2}\psi_{\text{o}})+\beta\frac{\partial\psi_{\text{o}}}{\partial
x}=-r\nabla^{2}\psi_{\text{o}}+\frac{C}{\rho
h}\nabla^{2}(\psi^{3}_{\text{a}}-\psi_{\text{o}}).$ (13)

The prognostic equations for the atmospheric and oceanic temperature fields,
using an energy balance scheme as in Barsugli and Battisti (1998), are

$\displaystyle\gamma_{\text{a}}\left(\frac{\partial T_{\text{a}}}{\partial
t}+J(\psi_{\text{a}},T_{\text{a}})-\sigma\omega\frac{p}{R}\right)$
$\displaystyle=-\lambda(T_{\text{a}}-T_{\text{o}})+\epsilon_{\text{a}}\sigma_{\text{B}}T_{\text{o}}^{4}-2\epsilon_{\text{a}}\sigma_{\text{B}}T_{\text{a}}^{4}+R_{\text{a}},$
(14) $\displaystyle\gamma_{\text{o}}\left(\frac{\partial
T_{\text{o}}}{\partial t}+J(\psi_{\text{o}},T_{\text{o}})\right)$
$\displaystyle=-\lambda(T_{\text{o}}-T_{\text{a}})-\sigma_{\text{B}}T_{\text{o}}^{4}+\epsilon_{\text{a}}\sigma_{\text{B}}T_{\text{a}}^{4}+R_{\text{o}}.$
(15)

The quartic terms in these equations are linearised by decomposing the
temperature fields around a spatially and temporally constant equilibrium
temperature, $T_{\text{a}}=T_{\text{a}}^{0}+\delta T_{\text{a}}$ and
$T_{\text{o}}=T_{\text{o}}^{0}+\delta T_{\text{o}}$, 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 $\delta T_{\text{a}}$ to the baroclinic streamfunction
$\theta_{\text{a}}\equiv(\psi_{\text{a}}^{1}-\psi_{\text{a}}^{3})/2$, more
specifically $\delta T_{\text{a}}=2f_{0}\theta_{\text{a}}/R$. Hence, the
remaining independent dynamical fields are the barotropic atmospheric
streamfunction field $\psi_{\text{a}}$, defined as
$\psi_{\text{a}}\equiv(\psi_{\text{a}}^{1}+\psi_{\text{a}}^{3})/2$, the
oceanic streamfunction field $\psi_{\text{o}}$, and the temperature anomalies
$\delta T_{\text{a}}$ and $\delta T_{o}$ of the atmosphere and the ocean. The
other parameters and variables that feature in the MAOOAM 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 fourth-order Runge-Kutta 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 $N_{x}$ and $N_{y}$ to specify the model resolution. If the
resolution of the ocean and the atmosphere are the same, the model resolution
is referred to as $N_{x}$x$N_{y}$; otherwise, it is denoted as atm.
$N_{x,\text{a}}$x$N_{y,\text{a}}$, oc. $N_{x,\text{o}}$x$N_{y,\text{o}}$.

Table 2: Variables and parameters in the MAOOAM equations Variable (units) | Description
---|---
$\psi_{\text{o}}$, $\psi_{\text{a}}$ (m2 s-1) | streamfunction of the ocean, atmosphere
$\omega=dp/dt$ (Pa s-1) | vertical velocity in pressure coordinates
$T_{\text{o}}$, $T_{\text{a}}$ (K) | temperature of the ocean, atmosphere; $T_{x}=T_{x}^{0}+\delta T_{x}$
$\delta T_{\text{o}}$, $\delta T_{\text{a}}$ (K) | temperature anomaly of the ocean, atmosphere
Parameter (units) | Description
$n=2L_{y}/L_{x}$ | meridional to zonal aspect ratio
$L_{y}=\pi L$ (km) | meridional extent of the domain
$f_{0}$ (s-1) | Coriolis parameter at 45∘ latitude
$\lambda$ (W m-2 K-1) | heat transfer coefficient at the ocean-atmosphere interface
$r$ (s-1) | friction coefficient at the bottom of the ocean layer
$C_{\text{o}}$, $C_{\text{a}}$ (W m-2) | insolation coefficient of the ocean, atmosphere
$k_{d}$ (s-1) | friction coefficient at ocean-atmosphere interface
$k_{d}^{\prime}$ (s-1) | friction coefficient between the atmospheric layers
$h$ (m) | depth of the ocean layer
$d=C/(\rho h)$ (s-1) | mechanical ocean-atmosphere coupling coefficient
$R$ (J kg-1 K-1) | gas constant of dry air
$L_{\text{R}}$ (km) | reduced Rossby deformation radius of the ocean
$\rho$ (kg m-3) | density of the ocean
$\sigma_{\text{B}}$ (W m2 K-4) | Stefan-Boltzmann constant
$\sigma$ (m2 s-2 Pa-2) | static stability of the atmosphere
$\beta$ (m-1 s-1) | Rossby parameter $\frac{df}{dy}$
$\gamma_{\text{o}}$, $\gamma_{\text{a}}$ (J m-2 K-1) | Specific heat capacity of the ocean layer, atmosphere
$T_{\text{o}}^{0}$, $T_{\text{a}}^{0}$ (K) | constant solution for the temperature of the ocean, atmosphere
$\epsilon_{\text{a}}$ | grey-body 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 as a dynamical system:

$\frac{d\vec{x}}{dt}={\vec{f}}(\vec{x},\alpha)$ (16)

where $\vec{x}$ is a vector containing the entire set of relevant variables
$\vec{x}$ = $(x_{1},...,x_{N})$ such as temperature or wind velocity,
projected on a relevant set of modes as described in Section 1. The function
${\vec{f}}$ is a nonlinear function of the variables $\vec{x}$ and $\alpha$
represents a set of parameters.

Let us consider a small perturbation along the trajectory, $\vec{x}(t)$,
generated by model (16), denoted $\delta\vec{x}(t)$. Provided that this
perturbation is sufficiently small (ideally infinitely small), its dynamics
can be described by the linearised equation,

$\frac{d\delta\vec{x}}{dt}={\left.\frac{\partial\vec{f}}{\partial\vec{x}}\right|}_{\vec{x}(t)}\delta\vec{x}$
(17)

and a formal solution can be written as

$\delta\vec{x}(t)={\bf M}(t,t_{0})\delta\vec{x}(t_{0})$ (18)

where the matrix $\bf{M}$ is referred to as the resolvent matrix or
propagator. This matrix $\bf{M}$ is responsible for the amplification or
contraction of the errors during the time period $t-t_{0}$. In order to get
information independent of the initial or final time, a limit
$(t-t_{0})\rightarrow\infty$ should be taken. Oseledets (1968, 2008)
demonstrates that provided that the system is ergodic, the following limit
exists for almost all initial conditions $\vec{x}(t_{0})=\vec{x}_{0}$,

$\lim_{t\rightarrow\infty}({\bf MM}^{\ast})^{1/2(t-t_{0})}=\Lambda_{x_{0}},$
(19)

where ${\bf M}^{\ast}$ is the adjoint of ${\bf M}$. The backward Lyapunov
exponents (Ershov and Potapov, 1998; Pazó et al., 2010) are then defined as
the natural logarithm of the eigenvalues of $\Lambda_{x_{0}}$. 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 ocean-atmosphere (Vannitsem and Lucarini, 2016).

If one or more LEs are positive, small errors on the initial conditions of the
system grow exponentially and the system is chaotic. In that case, the time
horizon of the system’s predictability is proportional to the inverse of the
largest Lyapunov exponent, $\frac{1}{\lambda_{1}}$. As this predictability
horizon is expressed in days for operational forecasting, we also express the
exponents $\lambda_{i}$ in units day-1. To translate the spectrum of LEs into
spatial scales of the instabilities in an unambiguous way, the CLVs must also
be determined. However, if there is scale-dependent dissipation, the largest
negative LEs are most likely to be associated with the smallest, most
dissipative scales.

The computation of the backward Lyapunov exponents follows the standard
algorithm of (Shimada and Nagashima, 1979; Benettin et al., 1980) based on the
Gram-Schmidt orthogonalisation.

  1. 1.

An ensemble ${\bf E}$ of $N$ perturbation vectors is randomly initialised.

  2. 2.

At every time step $t_{i}$, a matrix ${\bf P}_{i}$ that represents the linear
propagator from $t_{i-1}$ to $t_{i}$ is computed using the tangent linear
model along the model state trajectory. ${\bf P}_{i}$ is the equivalent of the
matrix $\bf{M}$ for a finite time difference $t_{i}-t_{i-1}$. We take into
account the numerical integration scheme when computing ${\bf P}_{i}$, by
evaluating the model Jacobian at all intermediate points of the scheme. We
have implemented the second- and fourth-order Runge-Kutta schemes, which
require two and four evaluations of the Jacobian per time step, respectively.

  3. 3.

As the model is integrated forward from time $t_{i}$ to $t_{i+b}$, the
corresponding linear propagator ${\bf P}_{i,i+b}$ is approximated by
multiplying the $b$ matrices, ${\bf P}_{i,i+b}={\bf P}_{i+b}\ldots{\bf
P}_{i+1}$. In the experiments that follow, we have chosen $b=1$.

  4. 4.

Every $b$ time steps, ${\bf E}$ is evolved with ${\bf P}_{i,i+b}$, and Gram-
Schmidt orthogonalised (using a $QR$-decomposition). The local Lyapunov
spectrum is computed from the diagonal of ${\bf R}$.

  5. 5.

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 Kaplan-Yorke or
Lyapunov dimension $D_{KY}$, which is an estimation of the information
dimension $D_{1}$. $D_{1}$ is known to be less than or equal to the capacity
or box-counting dimension $D_{0}$, also referred to as the fractal dimension
(Frederickson et al., 1983). $D_{KY}$ is defined as (Kaplan and Yorke, 1979):

$\displaystyle D_{KY}$
$\displaystyle=k+\frac{\lambda_{1}+\lambda_{2}+\ldots+\lambda_{k}}{|\lambda_{k+1}|},$
(20)

where $k$ is the highest index for which the sum of the largest $k$ Lyapunov
exponents is still strictly positive.

The second important property of the attractor is the Kolmogorov-Sinai or
metric entropy $h_{KS}$, 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,

$\displaystyle h_{KS}$ $\displaystyle\leq\sum_{\lambda_{i}>0}\lambda_{i},$
(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. FTLEs (Fujisaka, 1983; Abarbanel et
al., 1991) have been proposed to address such shortcomings, with the caveat
that they do not enjoy the extremely beneficial mathematical properties
(especially, norm-independence) 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 $\sigma^{j}_{m}$ over a time $m$ of
one backward Lyapunov exponent $\lambda_{j}$ and its statistics. As discussed
in (Schalge et al., 2013; Pazó et al., 2013; Laffargue et al., 2013), some
dynamical systems have the property that for a finite, but large $m$, the
fluctuations of their FTLEs can be described by a large deviation law (Kifer,
1990; Touchette, 2009). This is the case for Axiom A systems, and invoking the
chaotic hypothesis, extends to certain types of non Axiom A systems. For the
finite-time backward LEs and for a large $m$, we will verify the following
relation for the distribution $\mathcal{P}$ of the averages:

$\mathcal{P}\left(\sigma^{j}_{m}=x\right)\propto\exp\left(-m\,I_{j}(x)\right).$
(22)

$I_{j}(x)$ is the rate function, which is independent of $m$. The rate
function can be computed directly from this relation as

$I_{j}(x)=\lim_{m\rightarrow\infty}-\frac{1}{m}\,\log\left(\mathcal{P}\left(\sigma^{j}_{m}=x\right)\right).$
(23)

If $x$ represents a time series, we have to take the autocorrelation into
account. The FTLEs for the models under consideration have a non-zero
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 $T_{\text{decorr}}$. The time
$T_{\text{decorr}}$ is chosen to be the time lag when the autocorrelation
drops below $1/e$.

###  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 grid-point 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 spin-up 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 equator-to-pole temperature gradient $T_{EP}$
of 50 K and 60 K, respectively (with $\Delta T_{NS}=0$).

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 ocean-atmosphere coupling coefficient $d$ and the eddy
viscosity coefficients $\nu_{\text{a}}$ and $\nu_{\text{o}}$, as well as a
range of model resolutions.

Table 3: Model parameter values that are identical across all MAOOAM configurations used in this study. Parameter (unit) | Value | Parameter (unit) | Value
---|---|---|---
$n=2L_{y}/L_{x}$ | $1.5$ |  $L_{\text{R}}$ (km) | $19.93$
$L_{y}=\pi L$ (km) | $5.0\times 10^{3}$ |  $\rho$ (kg m-3) | $1000$
$f_{0}$ (s-1) | $1.032\times 10^{-3}$ |  $\sigma_{\text{B}}$ (W m2 K-4) | $5.6\times 5.610^{-8}$
$\lambda$ (W m-2 K-1) | $15.06$ |  $\sigma$ (m2 s-2 Pa-2) | $2.16\times 10^{-6}$
$r$ (s-1) | $1.0\times 10^{-7}$ |  $\beta$ (m-1 s-1) | $1.62\times 10^{-11}$
$C_{\text{o}}$ (W m-2) | $310$ |  $\gamma_{\text{o}}$ (J m-2 K-1) | $5.46\times 10^{8}$
$C_{\text{a}}$ (W m-2) | $C_{\text{o}}/3$ |  $\gamma_{\text{a}}$ (J m-2 K-1) | $1.0\times 10^{7}$
$k_{d}$ (s-1) | $3.0\times 10^{-6}$ |  $T_{\text{a}}^{0}$ (K) | $289.30$
$k_{d}^{\prime}$ (s-1) | $3.0\times 10^{-6}$ |  $T_{\text{o}}^{0}$ (K) | $301.46$
$h$ (m) | $136.5$ | $\epsilon_{\text{a}}$ | $0.7$
$R$ (J kg-1 K-1) | $287$ |  |

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 $10^{8}$ 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.
$b=1$. 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

This experiment corresponds to the setup of De Cruz et al. (2016). In addition
to the variables listed in Table 3, the mechanical ocean-atmosphere coupling
parameter $d$ is set to $1.1\times 10^{-7}$ s-1.

  * •

nodissip-reducedstress

For this “reduced-stress” experiment, the coupling parameter $d$ is reduced to
$4\times 10^{-8}\;\text{s}^{-1}$.

  * •

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:

$\displaystyle D_{o}$ $\displaystyle=\nu_{o}\nabla^{4}\psi_{o},$ (24)
$\displaystyle D_{a}$ $\displaystyle=\nu_{a}\nabla^{4}\psi_{a}.$ (25)

We adopt the values for the parameters $\nu_{o}$ and $\nu_{a}$ from Van der
Avoird et al. (2002), where they are estimated to be

$\displaystyle\nu_{o}$ $\displaystyle=1.76\times
10^{4}\;\text{m}^{2}\;\text{s}^{-1},$ (26) $\displaystyle\nu_{a}$
$\displaystyle=1\times 10^{5}\;\text{m}^{2}\;\text{s}^{-1}.$ (27)

Furthermore, $d$ is set to $1.1\times 10^{-7}\;s^{-1}$.

  * •

dissipationx10

In this experiment, $d=1.1\times 10^{-7}\;s^{-1}$, but the parameters
$\nu_{o}$ and $\nu_{a}$ are set to a higher value:

$\displaystyle\nu_{o}$ $\displaystyle=1.76\times
10^{5}\;\text{m}^{2}\;\text{s}^{-1},$ (28) $\displaystyle\nu_{a}$
$\displaystyle=1\times 10^{6}\;\text{m}^{2}\;\text{s}^{-1}.$ (29)

  * •

dissipation-reducedstress

This experiment has the same parameters as the “dissipation” experiment,
except for the coupling parameter $d$ which is reduced to $4\times
10^{-8}\;\text{s}^{-1}$.

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 subgrid-scale 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 the context of MAOOAM (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 200 largest LEs of the two different Lyapunov spectra
obtained in our experiments with PUMA. The averages were computed from a time
series of 25 years of daily finite-time LEs. We can estimate the size of the
attractor and by that estimate the degrees of freedom inherent to the
attractor with the Kaplan-Yorke dimension $D_{KY}$, as described in Section
2.1. The number of positive exponents and $D_{KY}$ are shown in the caption of
Fig. 1. Our findings confirm earlier results using two-layer QG models that
suggested an increase of $D_{KY}$ and the number of positive Lyapunov
exponents for a higher meridional temperature gradient (Lucarini et al., 2007;
Schubert and Lucarini, 2015).

Figure 1: The 200 largest LEs of the Lyapunov spectra of PUMA for the two
different setups with $\Delta T_{EP}=50K$ and $60K$. For $\Delta T_{EP}=50$K,
the number of positive Lyapunov exponents is 61 and $D_{KY}=172.6$. For
$\Delta T_{EP}=60$K, the number of positive Lyapunov exponents is 68 and
$D_{KY}=187.0$.

There are two very small exponents since the model setup is zonally symmetric
which in the limit of continuum creates an additional zero mode (see Schubert
and Lucarini (2015) for details). Otherwise, there are not many near-zero LEs
in PUMA. There is continuity between the time scales that characterise the QG
dynamics on the one hand, and faster, smaller-scale motions on the other hand.
This shows that the usual assumption of a clear time-scale 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
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 non-blocked states
of the mid-latitudes atmosphere.

We have computed the blocking rate employing the well established Tibaldi-
Molteni Index (Tibaldi and Molteni, 1990). We indeed find a higher blocking
rate for 50 K ($\approx 0.5$%) than for 60 K ($\approx 0.25$%). 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, 2016).

Next, the existence of a large deviation law for the FTLEs is verified, as
described in Section 2.2. The decorrelation time $T_{\text{decorr}}$ is
usually between 1 or 3 days. Therefore the rate function is computed for
averages of length $t_{\textrm{ave}}=m\cdot 3$ days.

In Figs. 3, 3, 7, 7, 9, 9, 5 and 5 the results for the rate functions are
shown for the fastest growing LE 1 (Figs. 3 and 3), fast-decaying LE 150
(Figs. 5 and 5), the positive and negative near-zero LEs 59 and 64 for $\Delta
T_{EP}=50K$ (Figs. 7 and 7) and near-zero LEs 66 and 71 for $\Delta
T_{EP}=60K$ (Figs. 9 and 9). Since Eq. 22 is needed to compute the rate
function, a long time series is necessary to estimate the distribution $P$
reliably.

Figure 2: Distributions and rate functions of $\lambda_{1}$, the fastest-
growing instability in PUMA, for $\Delta T_{EP}=50K$. The top panel shows the
different distributions and their kernel-smoothing approximation of
$\sigma^{1}_{m}$ where $t_{\textrm{ave}}$ is the respective $3\ m$. The bottom
panel shows a comparison of the rate functions, with the minimum moved to
zero.

Figure 3: Distributions and rate functions of $\lambda_{1}$, the fastest-
growing instability in PUMA, for $\Delta T_{EP}=60K$. Panels as in Fig. 3.

Figure 4: Distributions and rate functions of $\lambda_{150}$, a strongly
decaying direction in PUMA, for $\Delta T_{EP}=50K$. Panels as in Fig. 3.

Figure 5: Distributions and rate functions of $\lambda_{150}$, a strongly
decaying direction in PUMA, for $\Delta T_{EP}=60K$. Panels as in Fig. 3.

Figure 6: Distributions and rate functions of $\lambda_{59}$, a near-zero,
growing instability in PUMA, for $\Delta T_{EP}=50K$. Panels as in Fig. 3.

Figure 7: Distributions and rate functions of $\lambda_{64}$, a near-zero,
decaying instability in PUMA, for $\Delta T_{EP}=50K$. Panels as in Fig. 3.

Figure 8: Distributions and rate functions of $\lambda_{66}$, a near-zero,
growing instability in PUMA, for $\Delta T_{EP}=60K$. Panels as in Fig. 3.

Figure 9: Distributions and rate functions of $\lambda_{71}$ a near-zero,
decaying instability in PUMA, for $\Delta T_{EP}=60K$. Panels as in Fig. 3.

Our intent is to make at least a qualitative assessment of the convergence
rate for $m\rightarrow\ \infty$. The top panels of these figures show the
approximation of the respective distributions obtained via kernel density
estimation (Scott, 1979) of the distribution of the block-averaged LEs. The
bottom panels show the rate function for different $t_{\textrm{ave}}$ derived
from Eq. 22. A detailed investigation of other metrics, such as the
convergence of the variance, is provided in Appendix A.

We make the following observations. The graphs suggest a convergence of the
rate function for all LEs. Also, the rate functions’ shape is approximately
parabolic and the estimates of the rate functions appear to converge to the
asymptotic with a comparable speed regardless of the value of the
corresponding LE.

We interpret these results to stem from the lack of a clear-cut time-scale
separation in a purely atmospheric model like PUMA. This is in opposition to
what was originally speculated in Schubert and Lucarini (2015), where a
primitive-equation model was expected to feature a time-scale separation
visible in the Lyapunov spectrum. Such a time-scale separation would have been
an a posteriori justification of the filtering by the QG approximation.

We have shown that in a primitive-equation model with a high-dimensional phase
space of $\approx 60000$ the dimension of the attractor is small (order of
200) in comparison. Nevertheless, the unstable subspace can still be regarded
as a high-dimensional subspace with dimension of order of 50. 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 time-
scale 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 $\lambda_{1}$, 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 $\lambda_{1}$ 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 lower-order 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 scale-dependent dissipation term is present. The decrease in
$\lambda_{1}$ also appears to suggest an enhanced predictability for models
which have a larger ocean-atmosphere coupling parameter $d$, 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 low-order
36-variable model that lies at the basis of MAOOAM. Two distinct mechanisms
were identified to drive the increase in predictability with increasing $d$.
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 ocean-atmosphere coupling
above a critical value induces a sudden jump in predictability, associated
with the development of a slow coupled ocean-atmosphere mode (Vannitsem et
al., 2015; Vannitsem, 2017).

Figure 10: The largest Lyapunov exponent $\lambda_{1}$ of MAOOAM as a function
of the model resolution for the different experiments: nodissip (red pluses),
nodissip-reducedstress (green circles), dissipation (blue upward-pointing
triangles) dissipationx10 (purple diamonds), dissipation-reducedstress (orange
downward-pointing triangles).

Figures 11 to 14 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
near-zero 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 near-zero LEs in PUMA. Being a purely atmospheric model, PUMA’s Lyapunov
spectrum does not exhibit the large time-scale separation present in MAOOAM.
Indeed, the spectrum of PUMA bears more resemblance to that of the QG two-
layer model of Schubert and Lucarini (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 higher-resolution models are highly dissipative, hence
increasing the number of strongly negative Lyapunov exponents. The additional
positive and near-zero 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 $L_{\textit{Rh}}=\sqrt{\frac{U}{\beta}}$, requiring oceanic wavenumbers
as high as of 40–50.

Figure 11: Lyapunov spectra of MAOOAM for the “nodissip” experiment, for model
configurations from atm. 2x4, oc. 2x4 (red full line) up to atm. 11x11, oc.
11x11 (pink dash-dot-dotted line). Lyapunov exponents are ranked in decreasing
order, and the index of the smallest positive Lyapunov exponent is indicated
with a downward-pointing arrow for each model configuration. Figure 12:
Lyapunov spectra of MAOOAM for the “nodissip-reducedstress” experiment.
Colours and arrows as in Fig. 11. Figure 13: Lyapunov spectra of MAOOAM for
the “dissipation” experiment. Colours and arrows as in Fig. 11. Figure 14:
Lyapunov spectra of MAOOAM for the “dissipationx10” experiment. Colours and
arrows as in Fig. 11.

Figure 15 plots the Kaplan-Yorke dimension $D_{KY}$ as a function of the model
resolution. This shows that $D_{KY}$ is the highest for the models that do not
include the scale-dependent dissipation process (nodissip). A reduction in the
ocean-atmosphere coupling $d$ appears to slightly increase $D_{KY}$ for most
model resolutions, both in the case with and without scale-dependent
dissipation. The tenfold increase in the dissipation parameters $\nu_{o}$ and
$\nu_{a}$ (dissipationx10) results in the lowest values for $D_{KY}$, 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 $D_{KY}$ by the
number of dimensions $N$, as shown in Fig. 16. This shows that while $D_{KY}$
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 15: Kaplan-Yorke or Lyapunov dimension $D_{KY}$ of MAOOAM as a function
of the resolution for the different experiments. Colours as in Fig. 10. Figure
16: Kaplan-Yorke or Lyapunov dimension $D_{KY}$ of MAOOAM divided by the total
number of dimensions $N$, as a function of the resolution for the different
experiments. Colours as in Fig. 10. Figure 17: Kolmogorov-Sinai entropy
$h_{KS}$ of MAOOAM as a function of the resolution for the different
experiments. Colours as in Fig. 10.

Figure 17 shows the Kolmogorov-Sinai entropy $h_{KS}$ versus model resolution,
for the different experiments. The trends for the “nodissip” and “nodissip-
reducedstress” experiments appear to suggest that $h_{KS}$ would increase
unboundedly for increasing model resolution if a parametrization for the
scale-dependent dissipation is absent. The experiments that take this process
into account, paint a more realistic picture, with $h_{KS}$ 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. Indeed, the quantity $D_{KY}/N$, which approximates the relative
fraction of the attractor’s dimension, increases for increasing ocean
resolution, but decreases for increasing atmosphere resolution, as illustrated
in Fig. 19. This result deserves an extensive investigation by looking at the
properties of the CLVs.

Figure 18: Lyapunov spectra of MAOOAM for the “dissipation” experiment, for
different model configurations starting from atm. 6x6, oc. 6x6 (black full
line). The resolutions of the ocean (dark to light blue lines) and the
atmosphere (red to orange lines) are modified separately. Lyapunov exponents
are ranked in decreasing order, and the index of the smallest positive
Lyapunov exponent is indicated with a downward-pointing arrow for each model
configuration. Figure 19: Kaplan-Yorke or Lyapunov dimension $D_{KY}$ of
MAOOAM divided by the total number of dimensions $N$, as a function of the
model configuration.

As a final analysis, we have given a preliminary look at whether large
deviation laws can be established for the long-term 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 modes are weakly correlated. This
would suggest that one can construct the rate functions defining the large
deviations laws. The rate functions are shown in Fig. 20a-c) for the 351st LE.
Their convergence properties are investigated in Appendix A.2, and indicate
that we have not yet converged to the central limit theorem even for these
strongly damped modes. Additionally, the lagged time correlation of the near-
zero 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. Lagged-time
correlations are such that the available time series are not sufficiently long
to reach the asymptotic limit. This is even the case in the nodissip
simulation scenario, despite the larger value of the 1st LE and faster decay
of correlations; compare Figs. 21a)-c). This suggests that when many unstable
modes are present, disentangling their long-term properties requires very long
integrations, possibly as a result of geometrical quasi-degeneracies 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 ocean-atmosphere 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
low-order version of MAOOAM (Vannitsem and Lucarini, 2016), large deviation
laws cannot be established for the near-zero and positive FTLEs in the “9x9”
configuration.

Figure 20: Estimate of the rate function describing the large deviation law of
the 351st FTLE for MAOOAM with no dissipation (a), reference value for the
dissipation (b), and enhanced dissipation by a factor of 10 (c).

Figure 21: Estimate of the rate function describing the large deviation law of
the first FTLE for MAOOAM with no dissipation (a), reference value for the
dissipation (b), and enhanced dissipation by a factor of 10 (c).

[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
primitive-equation atmospheric model and a coupled ocean-atmosphere model.
Both systems suggest that high-dimensional 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 time-scale
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 quasi-geostrophic atmosphere-ocean model, which, by
definition, features a large time-scale separation between ocean and
atmosphere, and lacks a satisfactory representation of mesoscale and sub-
mesoscale processes. PUMA is an atmospheric-only primitive-equation model,
which can represent the faster, smaller-scale 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 built-in 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 quasi-geostrophic
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, Kelvin-Helmholtz instability,
etc. Despite this, we find that for the lower meridional temperature gradient
($\Delta T_{EP}=50K$) the spectrum features few near-zero negative exponents.
Interestingly, this might be related to the presence of blocking, but specific
studies with more robust dynamics between blocking and non-blocking situations
are necessary to clarify that. Additionally, the results suggest that the
FTLEs obey large deviation laws defining the predictability properties at long
time scales, including the near-zero 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 near-zero 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. Surprisingly, it is hard to find convergence for the rate functions
of the FTLEs, even for those 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 high-dimensional models to tackle various open problems. CLVs
allow to associate growth and decay rates to time-dependent 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 long-term 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 so-called
inertial manifold. The inertial manifold contains the effective finite-
dimensional dynamics of the system, which, we remind, is originally infinite-
dimensional 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 atmosphere-ocean 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 wave-
wave interactions (Speranza and Malguzzi, 1988) or regime switching like
blocking (Pelly and Hoskins, 2003). Additionally, CLVs will allow us to better
understand coupled ocean-atmospheric 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 quasi-geostrophic (QG) two-layer 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.uni-hamburg.de/en/arbeitsgruppen/theoretische-
meteorologie/modelle/sources/plasim.tgz. The source code for the latest
version of MAOOAM is available at http://github.com/Climdyn/MAOOAM. The
version of MAOOAM that was used to compute the Lyapunov exponents is archived
at https://doi.org/10.5281/zenodo.1198650.

The Lyapunov spectra of the different PUMA and MAOOAM configurations, that
were computed using the Benettin algorithm, are available as supplementary
material.

##  Appendix A Study of the convergence of the rate function

We have performed an additional analysis, that focuses on the convergence of
the standard deviation $\sigma$ of the FTLE as a function of the length
$t_{\textrm{ave}}$ used to compute the block averages. Indeed, $\sigma$ is
expected to scale according to $t_{\textrm{ave}}^{-\frac{1}{2}}$. Furthermore,
one expects the value of $\sigma^{2}\cdot t_{\textrm{ave}}$ to level off at
the value of the diffusion coefficient $D$ if there is convergence to the
central limit theorem. The diffusion coefficient $D$ is the inverse of the
second derivative of the rate function at its minimum, see e.g. (Kuptsov and
Politi, 2011).

###  A.1 PUMA

The scaling of $\sigma$ versus $t_{\textrm{ave}}$ has the expected behaviour
for both experiments conducted with PUMA, as shown in Fig. 22. This is also
apparent in the convergence of $\sigma^{2}\cdot t_{\textrm{ave}}$, shown in
Fig. 23. While the value of $D$ fluctuates, it has the right order of
magnitude.

Figure 22: Standard deviation $\sigma$ as a function of the block averaging
length $t_{\textrm{ave}}$ for different Lyapunov exponents of the PUMA model.
The Lyapunov index is shown in the legend. Panel (a) shows the results for a
temperature gradient $\Delta T_{EP}$ of 50K, panel (b) for 60K. The black
dashed line corresponds to $t_{\textrm{ave}}^{-\frac{1}{2}}$ scaling.

Figure 23: The metric $\sigma^{2}\cdot t_{\textrm{ave}}$ versus the diffusion
coefficient $D$, derived from the inverse of the second derivative at the
minimum of the rate function, as a function of the block averaging length
$t_{\textrm{ave}}$ for different Lyapunov exponents of the PUMA model. The
Lyapunov index is shown in the title. The top panels show the results for a
temperature gradient of 50K, the bottom panels for 60K.

###  A.2 MAOOAM

The decorrelation time of the near-zero FTLEs in MAOOAM is extremely long.
Accordingly, the time intervals $t_{\textrm{ave}}$, used to determine the rate
functions associated to the FTLE, are expected to be insufficient to robustly
define large deviations laws describing the statistics of the FTLEs. This is
confirmed by the behaviour of $\sigma$ versus $t_{\textrm{ave}}$ for the “9x9”
model configuration, shown in Fig. 24. A discrepancy between $D$ and
$\sigma^{2}\cdot t_{\textrm{ave}}$ is apparent for LE 100 in all experiments
shown in Fig. 25. Even for positive or strongly negative LEs, we are not yet
in the regime of convergence to the central limit theorem. Note however that
an integration time of 614 years was used to compute the Lyapunov spectra,
longer than the time series used in this analysis.

Figure 24: Standard deviation $\sigma$ as a function of the block averaging
length $t_{\textrm{ave}}$ for different Lyapunov exponents of the “9x9”
configuration of MAOOAM. The Lyapunov index is shown in the legend. Panel (a)
shows the results for the experiment without scale-dependent dissipation,
panel (b) corresponds to the reference value for dissipation, and panel (c)
shows the enhanced dissipation results. The black dashed line corresponds to
$t_{\textrm{ave}}^{-\frac{1}{2}}$ scaling.

Figure 25: The metric $\sigma^{2}\cdot t_{\textrm{ave}}$ versus the diffusion
coefficient $D$, derived from the inverse of the second derivative at the
minimum of the rate function, as a function of the block averaging length
$t_{\textrm{ave}}$ for different Lyapunov exponents of the “9x9” configuration
of MAOOAM. Note the logarithmic scale of the $y$-axis, in contrast to Fig. 23.
The Lyapunov index is shown in the title. The top panels show the results for
the experiment without scale-dependent dissipation, the centre panels
correspond to the reference value for dissipation, and the lower panels show
the results for an enhanced dissipation coefficient.

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.

The authors would like to thank S.G. Penny and an anonymous reviewer for
improving the manuscript with their constructive comments. 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

  * Abarbanel et al. (1991) Abarbanel, H. D., Brown, R., and Kennel, M. B.: Variation of Lyapunov exponents on a strange attractor, Journal of Nonlinear Science, 1, 175–199, 10.1007/BF01209065, 1991.
  * Barsugli and Battisti (1998) Barsugli, J. J. and Battisti, D. S.: The Basic Effects of Atmosphere-Ocean Thermal Coupling on Midlatitude Variability*, Journal of the Atmospheric Sciences, 55, 477–493, 10.1175/1520-0469(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, 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, 10.1016/S0370-1573(01)00025-4, 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, 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, 10.5194/npg-15-503-2008, 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, 10.1088/1751-8113/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, 10.1142/7351, 2010.
  * Charney (1947) Charney, J. G.: The Dynamics Of Long Waves In A Baroclinic Current, Journal of Meteorology, 4, 136–162, 10.1175/1520-0469(1947)004<0136:TDOLWI>2.0.CO;2, 1947.
  * Charney and Straus (1980) Charney, J. G. and Straus, D. M.: Form-Drag Instability, Multiple Equilibria and Propagating Planetary Waves in Baroclinic, Orographically Forced, Planetary Wave Systems, Journal of the Atmospheric Sciences, 37, 1157–1176, 10.1175/1520-0469(1980)037<1157:FDIMEA>2.0.CO;2, 1980.
  * De Cruz et al. (2016) De Cruz, L., Demaeyer, J., and Vannitsem, S.: The Modular Arbitrary-Order Ocean-Atmosphere Model: MAOOAM v1.0, Geoscientific Model Development, 9, 2793–2808, 10.5194/gmd-9-2793-2016, 2016.
  * Demaeyer and Vannitsem (2016) Demaeyer, J. and Vannitsem, S.: Stochastic parametrization of subgrid-scale processes in coupled ocean–atmosphere systems: benefits and limitations of response theory, Quarterly Journal of the Royal Meteorological Society, 143, 881–896, 10.1002/qj.2973, 2016.
  * Eady (1949) Eady, E. T.: Long Waves and Cyclone Waves, Tellus, 1, 33–52, 10.1111/j.2153-3490.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, 10.1103/RevModPhys.57.617, 1985\.
  * Ershov and Potapov (1998) Ershov, S. V. and Potapov, A. B.: On the concept of stationary Lyapunov basis, Physica D: Nonlinear Phenomena, 118, 167–198, 10.1016/S0167-2789(98)00013-X, 1998.
  * Fraedrich and Kirk (2005) Fraedrich, K. and Kirk, E.: The portable university model of the atmosphere (PUMA): Storm track dynamics and low-frequency variability, Meteorologische Zeitschrift, 14, 735–745, 10.1127/0941-2948/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, 10.1016/0022-0396(83)90011-6, 1983.
  * Fujisaka (1983) Fujisaka, H.: Statistical dynamics generated by fluctuations of local Lyapunov exponents, Progress of theoretical physics, 70, 1264–1275, 10.1143/PTP.70.1264, 1983.
  * Gallavotti (2014) Gallavotti, G.: Nonequilibrium and Irreversibility, Theoretical and Mathematical Physics, Springer, New York, 10.1007/978-3-319-06758-2, 2014\.
  * Gallavotti and Cohen (1995) Gallavotti, G. and Cohen, E. G. D.: Dynamical Ensembles in Nonequilibrium Statistical Mechanics, Physical Review Letters, 74, 2694–2697, 10.1103/PhysRevLett.74.2694, 1995.
  * Gallavotti and Lucarini (2014) Gallavotti, G. and Lucarini, V.: Equivalence of Non-equilibrium Ensembles and Representation of Friction in Turbulent Flows: The Lorenz 96 Model, Journal of Statistical Physics, 156, 1027–1065, 10.1007/s10955-014-1051-6, 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, 10.1016/0375-9601(91)90012-W, 1991.
  * Giering and Kaminski (2003) Giering, R. and Kaminski, T.: Applying TAF to generate efficient derivative code of Fortran 77-95 programs, PAMM, 2, 54–57, 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, 10.1103/PhysRevLett.99.130601, URL http://arxiv.org/abs/0706.0510, 2007.
  * 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, 10.1175/1520-0477(1994)075<1825:APFTIO>2.0.CO;2, 1994.
  * Holden et al. (2014) Holden, P. B., Edwards, N. R., Garthwaite, P. H., Fraedrich, K., Lunkeit, F., Kirk, E., Labriet, M., Kanudia, A., and Babonneau, F.: PLASIM-ENTSem v1.0: a spatio-temporal emulator of future climate change for impacts assessment, Geoscientific Model Development, 7, 433–451, 10.5194/gmd-7-433-2014, 2014\.
  * Holton (2004) Holton, J. R.: An Introduction To Dynamic Meteorology, Academic Press, San Diego, California, USA, fourth edn., 2004.
  * Kalnay (2003) Kalnay, E.: Atmospheric Modeling, Data Assimilation, and Predictability, Cambridge University Press, New York, USA, 2003.
  * Kaplan and Yorke (1979) Kaplan, J. L. and Yorke, J. A.: Chaotic behavior of multidimensional difference equations, in: Functional Differential Equations and Approximation of Fixed Points: Proceedings, Bonn, July 1978, edited by Peitgen, H.-O. and Walther, H.-O., pp. 204–227, Springer-Verlag, Berlin, 1979.
  * Kifer (1990) Kifer, Y.: Large deviations in dynamical systems and stochastic processes, Transactions of the American Mathematical Society, 321, 505–524, 10.1090/S0002-9947-1990-1025756-7, 1990.
  * Klein (2010) Klein, R.: Scale-Dependent Models for Atmospheric Flows, Annual Review of Fluid Mechanics, 42, 249–274, 10.1146/annurev-fluid-121108-145537, 2010.
  * Kuptsov and Politi (2011) Kuptsov, P. V. and Politi, A.: Large-deviation approach to space-time chaos, Physical Review Letters, 107, 114 101, 10.1103/PhysRevLett.107.114101, 2011\.
  * 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, 10.1088/1751-8113/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, 10.1175/1520-0469(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, 10.1007/978-94-009-0533-7_1, 1990.
  * Lorenz (1963) Lorenz, E. N.: Deterministic Nonperiodic Flow, Journal of the Atmospheric Sciences, 20, 130–141, 10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2, 1963\.
  * Lucarini and Fraedrich (2009) Lucarini, V. and Fraedrich, K.: Symmetry breaking, mixing, instability, and low-frequency variability in a minimal Lorenz-like system, Physical Review E, 80, 026 313, 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, 10.1029/2009RG000323, 2011.
  * Lucarini et al. (2007) Lucarini, V., Speranza, A., and Vitolo, R.: Parametric smoothness and self-scaling of the statistical properties of a minimal climate model: What beyond the mean field theories?, Physica D: Nonlinear Phenomena, 234, 105–123, 10.1016/j.physd.2007.07.006, 2007.
  * Lucarini et al. (2010) Lucarini, V., Fraedrich, K., and Lunkeit, F.: Thermodynamics of climate change: generalized sensitivities, Atmospheric Chemistry and Physics, 10, 9729–9737, 10.5194/acp-10-9729-2010, URL http://www.atmos-chem-phys.net/10/9729/2010/, 2010.
  * 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, 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, 10.1007/s10955-016-1506-z, 2017.
  * Manneville (1985) Manneville, P.: Liapounov exponents for the Kuramoto-Sivashinsky model, in: Macroscopic Modelling of Turbulent Flows, vol. 230, pp. 319–326, Springer Berlin Heidelberg, Berlin, Heidelberg, 10.1007/3-540-15644-5_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, 10.1007/3-540-60188-0_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, 10.1175/JPO2770.1, 2005.
  * Nese and Dutton (1993) Nese, J. M. and Dutton, J. A.: Quantifying predictability variations in a low-order occan-atmosphere model: a dynamical systems approach, Journal of Climate, 6, 185–204, 10.1175/1520-0442(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, 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, 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, 10.1175/2008JAS2781.1, 2009.
  * Oseledets (2008) Oseledets, V.: Oseledets theorem, Scholarpedia, 3, 1846, 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, 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, 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.: Spatio-temporal evolution of perturbations in ensembles initialized by bred, Lyapunov and singular vectors, Tellus, Series A: Dynamic Meteorology and Oceanography, 62, 10–23, 10.1111/j.1600-0870.2009.00419.x, 2010.
  * Pazó et al. (2013) Pazó, D., López, J. M., and Politi, A.: Universal scaling of Lyapunov-exponent fluctuations in space-time chaos, Physical Review E, 87, 062 909, 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, 10.1175/1520-0469(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, 10.1175/1520-0469(2003)060<0743:ANPOB>2.0.CO;2, 2003.
  * Pesin (2004) Pesin, Y.: Lectures on Partial Hyperbolicity and Stable Ergodicity, European Mathematical Society, Zurich, Switzerland, 2004.
  * Pierini (2011) Pierini, S.: Low-frequency variability, coherence resonance, and phase selection in a low-order model of the wind-driven ocean circulation, Journal of Physical Oceanography, 41, 1585–1604, 10.1175/JPO-D-10-05018.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, 10.1002/2012RG000419, 2014.
  * Ragone et al. (2016) Ragone, F., Lucarini, V., and Lunkeit, F.: A new framework for climate sensitivity and prediction: a modelling perspective, Climate Dynamics, 46, 1459–1471, 10.1007/s00382-015-2657-3, 2016.
  * Ruelle (1979) Ruelle, D.: Ergodic theory of differentiable dynamical systems, Publications mathématiques de l’IHÉS, 50, 27–58, 10.1007/BF02684768, 1979\.
  * Ruelle (2009) Ruelle, D.: A review of linear response theory for general differentiable dynamical systems, Nonlinearity, 22, 855–870, 10.1088/0951-7715/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, 87, 10.1103/PhysRevE.87.052113, URL https://mi-pub.cen.uni-hamburg.de/fileadmin/files/forschung/theomet/docs/pdf_2013/2013-Schalge-fluct-theor-PRE.pdf, 2013\.
  * Schneider (2006) Schneider, T.: The General Circulation of the Atmosphere, Annual Review of Earth and Planetary Sciences, 34, 655–688, 10.1146/annurev.earth.34.031405.125144, 2006.
  * Schubert and Lucarini (2015) Schubert, S. and Lucarini, V.: Covariant Lyapunov vectors of a quasi-geostrophic baroclinic model: Analysis of instabilities and feedbacks, Quarterly Journal of the Royal Meteorological Society, 141, 3040–3055, 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, 10.1002/qj.2808, 2016\.
  * Scott (1979) Scott, D. W.: On optimal and data-based histograms, Biometrika, 66, 605–610, 1979\.
  * 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, 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: Quasi-geostrophic, Two-Layer Model Atmosphere, Journal of the Atmospheric Sciences, 45, 3046–3062, 10.1175/1520-0469(1988)045<3046:TSPOAZ>2.0.CO;2, 1988.
  * Sprott (2010) Sprott, J. C.: Elegant Chaos, World Scientific, 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, 10.3402/tellusa.v9i3.9111, 1957.
  * Tibaldi and Molteni (1990) Tibaldi, S. and Molteni, F.: On the operational predictability of blocking, Tellus A, 42, 343–365, 10.1034/j.1600-0870.1990.t01-2-00003.x, 1990.
  * Touchette (2009) Touchette, H.: The large deviation approach to statistical mechanics, Physics Reports, 478, 1–69, 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, 10.1175/1520-0469(1998)055<0390:POLVAS>2.0.CO;2, 1998.
  * Trevisan et al. (2010) Trevisan, A., D’Isidoro, M., and Talagrand, O.: Four-dimensional 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 large-scale circulation, Cambridge University Press, New York, USA, 2006.
  * Van der Avoird et al. (2002) Van der Avoird, E., Dijkstra, H., Nauw, J., and Schuurmans, C.: Nonlinearly induced low-frequency variability in a midlatitude coupled ocean-atmosphere model of intermediate complexity, Climate Dynamics, 19, 303–320, 10.1007/s00382-001-0220-x, 2002.
  * Vannitsem (2017) Vannitsem, S.: Predictability of large-scale atmospheric motions: Lyapunov exponents and error dynamics, Chaos: An Interdisciplinary Journal of Nonlinear Science, 27, 32 101, 10.1063/1.4979042, 2017.
  * Vannitsem and De Cruz (2014) Vannitsem, S. and De Cruz, L.: A 24-variable low-order coupled ocean–atmosphere model: OA-QG-WS v2, Geoscientific Model Development, 7, 649–662, 10.5194/gmd-7-649-2014, 2014.
  * Vannitsem and Lucarini (2016) Vannitsem, S. and Lucarini, V.: Statistical and Dynamical Properties of Covariant Lyapunov Vectors in a Coupled Atmosphere-Ocean Model - Multiscale Effects, Geometric Degeneracy, and Error Dynamics, Journal of Physics A: Mathematical and Theoretical, 49, 224 001, 10.1088/1751-8113/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, 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, 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, 10.1175/1520-0469(1997)054<0347:LVAEGP>2.0.CO;2, 1997.
  * Vannitsem et al. (2015) Vannitsem, S., Demaeyer, J., De Cruz, L., and Ghil, M.: Low-frequency variability and heat transport in a low-order nonlinear coupled ocean-atmosphere model, Physica D: Nonlinear Phenomena, 309, 71–85, 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, 10.1111/j.1600-0870.2007.00234.x, 2007.
  * Xu and Paul (2016) Xu, M. and Paul, M. R.: Covariant Lyapunov vectors of chaotic Rayleigh-Bénard convection, Physical Review E, 93, 062 208, 10.1103/PhysRevE.93.062208, 2016.
  * Yamada and Ohkitani (1988) Yamada, M. and Ohkitani, K.: Lyapunov spectrum of a model of two-dimensional turbulence, Physical Review Letters, 60, 983–986, 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, 10.1088/1751-8113/46/25/254015, 2013.