[1]Guilherme Tegoni Goedert

[1]School for Applied Mathematics, Getúlio Vargas Foundation, Praia de
Botafogo, 190, Rio de Janeiro, 22250-900, Rio de Janeiro, Brazil

2]Department of Physics and INFN, University of Rome, Tor Vergata, Via della
Ricerca Scientifica 1, Rome, 00133, Lazio, Italy

# Instanton-based Importance Sampling for Extreme Fluctuations in a Shell
Model for Turbulent Energy Cascade

<EMAIL_ADDRESS>Luca Biferale<EMAIL_ADDRESS>*  [

###### Abstract

Many out-of-equilibrium flows present non-Gaussian fluctuations in physically
relevant observables, such as energy dissipation rate. This implies extreme
fluctuations that, although rarely observed, have a significant phenomenology.
Recently, path integral methods for importance sampling have emerged from
formalism initially devised for quantum field theory and are being
successfully applied to the Burgers equation and other fluid models. We
proposed exploring the domain of application of these methods using a Shell
Model, a dynamical system for turbulent energy cascade which can be
numerically sampled for extreme events in an efficient manner and presents
many interesting properties. We start from a validation of the instanton-based
importance sampling methodology in the heat equation limit. We explored the
limits of the method as non-linearity grows stronger, finding good qualitative
results for small values of the leading non-linear coefficient. A worst
agreement between numerical simulations of the whole systems and instanton
results for estimation of the distribution’s flatness is observed when
increasing the nonlinear intensities.

###### keywords:

turbulence, shell models, extreme events, importance sampling, instantons

##  1 Introduction

Turbulent flow is a complex system where coherent structures are formed in
time and space and interact on a wide range of spatial scales, as energy is
introduced in large scales (between kilometers for atmospheric flows and
centimeters for a cup of coffee) but is dissipated only at scales many orders
of magnitude smaller. In between energy injection and dissipation, there is a
wide range of scales where these structures, eddies, exhibit self-similar
behavior and a well-recorded power-law spectrum for energy; this is called the
inertial range [1].

Of particular interest to the investigation of this phenomenology are shell
models for turbulent energy cascade. These models can present highly non-
trivial and intermittent behavior analogous to fully developed turbulence [1,
2, 3, 4, 5]. Nevertheless, these dynamical systems are numerically and
analytically more tractable than the Navier-Stokes equations, serving as a
very interesting class of models for testing ideas and methods in fluid
dynamics.

In studying out-of-equilibrium micro- and macroscopic fluids, strongly non-
Gaussian and intermittent fluctuations are frequently observed. This implies
that abrupt changes in fluid velocity occur more often and with larger effects
than would be observed if Gaussian statistics were followed. These are often
related to sparse but intense coherent structures both at large scales (e.g.
atmospheric fronts, tornadoes) and small scales (influencing dissipation of
energy in the atmosphere or the ocean [6, 7], or the rain formation process
[8]).

However, the intrinsic rarity of these extreme events makes their statistical
assessment a great challenge. To investigate extreme events in the energy
cascade of classical hydrodynamic turbulence, we propose the employment of
Instanton-based methods for the importance sampling of intense fluctuations
and explore this methodology in Shell Models, a class of dynamical systems
which mimic turbulent energy cascade while being an efficient playground for
sampling of extreme event numerically.

Instanton-based methods for importance sampling have been successfully tested
for the Burgers equations [9, 10, 11, 12]. Furthermore, instantons as
statistically dominant field configurations have been used to describe
circulation intermittency in Navier-Stokes [13, 14]. Here, we explore the
simpler but dynamically rich setting of shell models for hydrodynamics
turbulence, which have been used also as testing ground for modeling
magnetohydrodynamics [15, 16, 17], passive scalars [18, 19, 20] and convective
turbulence [21, 22], among others.

##  2 The Densniansky-Novikov shell model

Shell models are systems of (infinitely many) coupled ordinary differential
equations, modeled after the Navier-Stokes equations in Fourier space and
which preserve many of this system’s nonlinear properties and symmetries. One
can create shell models by discretizing the Fourier space into concentric
shells with radii arranged in a geometric sequence ($k_{n}=k_{0}h^{n}$), and
keeping only one variable $u_{n}(t)$ -or a few variables-representing each
shell. This is a further simplification with respect to projecting the
equation into a logarithmic lattice [23], and this greatly decreases the
number of degrees of freedom required to simulate solutions (to the order of
the logarithm of the Reynolds number, rather than a power of it).

The real-valued Desniansky-Novikov (D.N.) shell model, generalized to extend
the Obukhov model and with a stochastic forcing, reads [24, 5]

$\begin{split}\dfrac{du_{n}}{dt}=&\ G_{n}[u]+f_{n},\ \ \left\langle
f_{n}(t)f_{m}(t^{\prime})\right\rangle=\chi_{nm}\delta(t-t^{\prime})\\\
G_{n}[u]=&\ c_{1}(k_{n}u_{n-1}^{2}-k_{n+1}u_{n}u_{n+1})\\\ +&\
c_{2}(k_{n}u_{n-1}u_{n}-k_{n+1}u_{n+1}^{2})-\nu k_{n}^{2}u_{n},\end{split}$
(1)

where $c_{1}$, $c_{2}$ and $\nu$ are real constant parameters, with $\nu\geq
0$ corresponding to the viscosity and $\chi_{nm}$ is the correlation between
the forcing applied to different shells. Although this is an infinite-
dimensional dynamical system, in practice we need to truncate the shell space
to $n\in\left\\{0,1,...,N\right\\}$ while imposing boundary conditions
$u_{-1}=u_{N+1}=0$ and choosing $N$ large enough to include integral,
inertial, and viscous scales [1]. This model was constructed to conserve total
energy $E=\frac{1}{2}\sum_{n}u_{n}^{2}$ in inviscid and unforced settings.

Among the key similarities between Shell Models and fully developed
turbulence, is the existence of a turbulent energy cascade “à la Kolmogorov”.
An example of energy spectrum is shown in Figure 1 for the Desniansky-Novikov
shell model (1) where we verify that
$E(k_{n})=\langle|u_{n}|^{2}\rangle\propto k_{n}^{-2/3}$, which is the
equivalent of the $-5/3$ power-law spectrum in three dimensional isotropic
turbulence considering the extra factor included by the thickness of each
shell. Here and hereafter time average is denoted with $\langle\cdot\rangle$.
Moreover, shell models also present highly intermittent bursts in its energy
dissipation rate $\epsilon=\nu\sum_{n}k_{n}^{2}u_{n}^{2}$, with many episodes
an order of magnitude or more higher than the average value. Nonetheless,
these rare events are a lot easier to sample in Shell Models than in the
Navier-Stokes equation highlighting the usefulness of these simpler models,
specially in the study of extreme events.

Figure 1: Mean energy values per shell $E(k_{n})$ in the D.N. Shell Model (1).
Parameters $c_{1}=1$, $c_{2}=0$ regulate the energy transfer between scales,
while we vary the values of viscosity $\nu$.

The D.N. shell model was chosen for this study due to its simplicity as well
as its different asymptotic behaviors when changing the control parameters
$c_{1}$ and $c_{2}$ [5, 24]. In the absence of stochastic forcing, the
solutions of this shell model tend to fixed-point attractors compatible with
the power-law scaling found in the energy spectrum. As soon as a stochastic
forcing is introduced, the system is driven out of equilibrium while still
preserving the scaling predicted by the energy spectrum (see Figure 2). At the
same time, intermittent but large fluctuations appear in the high wavenumber
shells, translating into bursts in the energy dissipation rate much larger
than its average value which are analogous to those observed in turbulent
flows. Even in a weakly nonlinear setting (see later Section 4.3), the
distribution of the energy dissipation rate showed heavy tails as seen in
Figure 3.

Figure 2: Solution for the stochastically forced D.N. shell model (1) on the
left and corresponding time series for the energy dissipation rate
$\epsilon=\nu\sum k_{n}^{2}u_{n}^{2}$ on the right; parameters used:
$c_{1}=c_{2}=0.5$, $\nu=10^{-4}$ and $\chi_{nm}=k_{n}^{-6}\delta_{nm}$. Figure
3: PDF of energy dissipation rate $\epsilon$ for a large dataset with $10^{8}$
samples, in the weakly linear setting explored for validation of our
methodology ($c_{1}=10^{-3}$, $\nu=10^{-2}$, $c_{2}=0$ and
$\chi_{nm}=k_{n}^{-6}\delta_{nm}$). The average value
$\langle\epsilon\rangle=0.5166$,

##  3 Instanton calculation and importance sampling

Strong non-Gaussian and intermittent fluctuations are observed in many out-of-
equilibrium micro and macroscopic fluids, implying that abrupt changes in
fluid velocity happen more often and with larger effects than would be
observed if Gaussian statistics were followed. These large gradients are often
related to sparse but intense coherent structures at both large scales (e.g.,
atmospheric fronts, tornadoes) and small scales (influencing dissipation of
energy in the atmosphere or the ocean [6, 7], or the rain formation process
[8]).

Despite the importance of such phenomena, their rarity leads to exceeding
difficulties in their observation, in the real world or in simulations.
Therefore, it is key to have numerical tools capable of quantitatively
assessing the statistical weight of these coherent or quasicoherent spatio-
temporal instanton-type solutions on the statistical ensemble.

Only very recently, new numerical techniques based on constrained Monte Carlo
approaches emerging from the lattice QCD community have been proposed to
follow this strategy in simple one-dimensional stochastic partial differential
equations [25, 26, 10, 27, 28].

###  3.1 Path integral formulation of stochastic differential equations

These techniques are mainly based on path integral Martin-Siggia-Rose-Janssen-
de Dominicis (MSRJD) formalism [29, 30, 31]. The general set up can be
described by a differential equation for a field $u(x,t)$, forced by some
normally distributed random noise $f$ which is white in time and has some
desired correlation $\chi$,

$\begin{split}&\dot{u}=N[u]+f(t),\\\ &\left\langle
f(t),f(t^{\prime})\right\rangle=\chi\delta(t-t^{\prime}).\end{split}$ (2)

In the context of shell models (1), $u=[u_{0},u_{1},...]$,
$N[u]=[G_{0}[u],G_{1}[u],...]$ and the stochastic forcing $f$ is also defined
element-wise for the discrete shell space $k_{n}$ as $f(t)=[f_{0},f_{1},...]$
with $\left\langle
f_{n}(t)f_{m}(t^{\prime})\right\rangle=\chi_{n,m}\delta(t-t^{\prime})$. The
MSRJD formalism allows for the direct evaluation of the mean value of a
desired observable $\mathcal{O}[u]$ by using a path integral written as the
Onsager-Machlup functional [32, 33]:

$\begin{split}&\left\langle\mathcal{O}[u]\right\rangle=\int\mathcal{D}u\mathcal{O}[u]J[u]e^{-S_{\mathcal{L}}[u,\dot{u}]},\\\
&S_{\mathcal{L}}[u,\dot{u}]=\dfrac{1}{2}\int
dt\left\langle\dot{u}-N[u],\chi^{-1}(\dot{u}-N[u])\right\rangle,\end{split}$
(3)

where $J[u]$ is the Jacobian for the transformation between a noise
realisation and velocity $f\rightarrow u$, $\mathcal{D}u$ is the element over
which functional integration is to be performed and
$S_{\mathcal{L}}[u,\dot{u}]$ is the Lagrangian action defined by this
formalism for the system (LABEL:eq:dynamics).

We perform the Hubbard-Stratonovich transformation [34, 35] on (3), therefore,
defining new momenta variables $p(x,t)$ and circumventing the dependency on
the inverse correlator $\chi^{-1}$. The resulting path integral is known as
the MSRJD response functional

$\begin{split}&\left\langle\mathcal{O}[u]\right\rangle=\int\mathcal{D}u\mathcal{D}p\mathcal{O}[u]J[u]e^{-S_{\mathcal{H}}[u,p]},\\\
&S_{\mathcal{H}}[u,p]=\int dt\left[\left\langle
p,\dot{u}-N[u]\right\rangle-\dfrac{1}{2}\left\langle p,\chi
p\right\rangle\right].\end{split}$ (4)

In the general context of shell models, of particular interest is the use of
observable that constrains certain functionals $O[u(t=0)]$ of the fields to
desired values

$\mathcal{O}[u]=\delta(O[u(t=0)]-a).$ (5)

The time at which the observable (5) is evaluated is arbitrary when dealing
with autonomous systems and we adopt the convention where the observation time
is shifted to the origin and the system is integrated over negative values of
$t$. In this case, the system may be rewritten in terms of a constrained
action $S_{\lambda}[u,p]$

$\begin{split}&\left\langle\mathcal{O}[u]\right\rangle=\int\mathcal{D}u\mathcal{D}p\dfrac{1}{2\pi}\int_{-\infty}^{\infty}d\lambda
e^{-S_{\lambda}[u,p]},\\\
&S_{\lambda}[u,p]=\lambda\left(O[u]-a\right)+S_{\mathcal{H}}[u,p].\end{split}$
(6)

By redefining the problem from a stochastic differential system to a path
integral formulation, one paves the way to the use of completely different
methodological tools. Let us notice that Hybrid Monte Carlo (HMC) [36, 37, 38]
has also been successfully implemented via modification of the Onsager-Machlup
functional to define an artificial Hamiltonian whose evolution is used to find
candidates for states in a Metropolis–Hastings algorithm; this method is
explored in [10] for the stochastically forced Burgers equations.

###  3.2 Instanton equations for shell models

Instantons are field configurations that are extrema of the action, that is,
they are the largest contributors to the evaluation of the mean value of an
observable in the path integral formulation. Therefore, they can be found as
solutions to the equations obtained via a saddle-point approximation. Here, we
are interested to realizations where a given function attains desired value.
An instanton is therefore found as a solution to the system of equations
derived through a saddle-point approximation of the action in (6),

$\dfrac{\delta S_{\lambda}}{\delta u}=0,\quad\dfrac{\delta S_{\lambda}}{\delta
p}=0.$ (7)

for a system of the form (LABEL:eq:dynamics), the conditions (7) lead to a
coupled set of equations between the original field and the momenta:

$\begin{split}\dot{u}&=N[u,x]+\chi\ast p,\\\
\dot{p}&=-\left(\nabla_{u}N[u]\right)^{T}p+\lambda\nabla_{u}O[u(t)]\delta(t),\end{split}$
(8)

where $\ast$ stands for the convolution product and $\delta$ is a Dirac delta
that is nonzero at the observation time defined in (5). The idea is that we
want to observe an extreme fluctuation at this observation time, giving it
enough time to develop [12]. Therefore, this system is defined on $t\in(-T,0]$
where $T>0$ is chosen large enough to satisfy condition $u$, $p\rightarrow 0$
as $t\rightarrow-\infty$. We observe that due to the Dirac delta in the second
equation of (8), the term $\lambda\nabla_{u}O[u(t)]$ is effective only at the
observation time $t=0$ at the right extreme of the time interval; this term
can therefore be removed from the equation and considered as a boundary
condition. The MSRJ response functional for a generic shell model is:

$\begin{split}S_{\mathcal{H}}[u,p]=\int
dt\left(\sum_{m=0}^{N}p_{m}(\dot{u}_{m}-G_{m}[u])\right.\\\
\left.-\dfrac{1}{2}\sum_{m=0}^{N}\sum_{l=0}^{N}p_{m}\chi_{ml}p_{l}\right).\end{split}$
(9)

The instanton equations then become ($s$ being the number of shells between
the farthest interacting neighbouring shell in the model)

$\dot{u}_{n}=G_{n}[u]+\sum_{l=0}^{N}\chi_{nl}p_{l}\
,\quad\dot{p}_{n}+\sum_{m=n-s}^{n+s}p_{m}\dfrac{\partial G_{m}}{\partial
u_{n}}=0,$ (10)

defined on $t\in(-T,0]$ which leads to the boundary condition
$p_{n}(0)=\lambda\nabla_{u}O[u(0)]$. We choose $T>0$ large enough to satisfy
the condition $u$, $p\rightarrow 0$ on the left boundary. For the real-valued
Desniansky-Novikov model (1), the instanton equations become

$\begin{split}\dot{u}_{n}+\nu
k_{n}^{2}u_{n}&=c_{1}k_{n}\left(u_{n-1}^{2}-hu_{n}u_{n+1}\right)\\\
&+c_{2}k_{n}\left(u_{n-1}u_{n}-hu_{n+1}^{2}\right)+\sum_{l=0}^{N}\chi_{nl}p_{l},\\\
\dot{p}_{n}-\nu
k_{n}^{2}p_{n}&=-c_{1}(2k_{n+1}u_{n}p_{n+1}-k_{n+1}u_{n+1}p_{n}\\\
&-k_{n}u_{n-1}p_{n-1})-c_{2}(k_{n+1}u_{n+1}p_{n+1}\\\
&+k_{n}u_{n-1}p_{n}-2k_{n}u_{n}p_{n-1}),\end{split}$ (11)

with the extra constraints $p_{n}(0)=-i\lambda\nabla_{u}O[u(0)]$ and
$u,p\rightarrow 0$ as $t\rightarrow-\infty$.

###  3.3 Chernykh-Stepanov iterative method

Though simple in nature, instanton equations generally pose an interesting
numerical challenge: if $N[u]$ contains a linear dissipative term, as is
usually the case in viscous systems, the equations for $u$ and $p$ are
unstable if integrated together in the same direction of time (see the
different signs of the viscous terms in system (11)) [12, 33]. As such, we
need to employ a specialized iterative method that splits the dynamics of each
field and solves each along its stable direction in time, as described in
Figure 4 [39, 12, 33, 9, 40]. Note that the dynamics (11) of one field
includes a coupling term with the other; in order to evaluate this coupling,
we need to solve both equations at the same instant of time, and we use the
approximations from the last iteration of the method at the same instant to
evaluate the coupling terms during integration. Note that the term
accompanying the Dirac delta in the evolution of $p$ can be considered a
simple initial condition $p(0)=-i\lambda\nabla_{u}O[u]$. Observe that the
constrained value $a$ of the functional $O[u]$ does not take part in the
dynamics; the instanton configuration for a given dynamics and functional
depends only on the parameter $\lambda$ that arises from the Fourier transform
of the constraint (5) in order to incorporate it into the action (6). The
constraint value is only used in the calculation of the action $S_{\lambda}$,
which is used to calculate the PDF for the observable studied up to a
prefactor as $\exp(-S_{\lambda})$.

Figure 4: Graphic representation of the Chernykh-Stepanov iterative method.

##  4 Solving instanton equations in shell models

###  4.1 Locality of observables and forces

Observables and forces play a key role in the formulation of instanton
equations and even more so in their solution. It is the spatial correlator
$\chi$ of the forces that couple the dynamics of $u$ to the configuration of
the auxiliary variable $p$. Also, it is the u-gradient of the constraint
functional $O[u]$ that injects energy into the system by serving as the
initial condition for $p$ as it is solved backward in time. However, not all
choices of force and observable are compatible in order to develop non trivial
instanton solutions.

In order to initiate the Chernykh-Stepanov iterative scheme, we begin with the
fields at rest with $u(t)=p(t)=0$ everywhere except the right boundary where
$\lambda\nabla_{u}O[u(0)]$ is enforced for a fixed value of $\lambda$. The
system is excited by the right boundary condition and energy it introduces is
initially propagated due to the term $\sum_{l}\chi_{nl}p_{l}$. Let us suppose
that an observable local to shell $n_{o}$ is chosen (i.e., the observable is a
function only of the component $n_{o}$ of $u$), while a different shell
$n_{f}$ is the only one forced. Since $G_{n}[u]$ is composed by a linear
viscous term and a sum of bilinear terms of shifts of $u$, it is easy to see
that in the first iteration of the Chernykh-Stepanov method, only
$p_{n_{o}}\neq 0$. However, if a force is local and $n_{o}\neq n_{f}$, then
$\dot{u}_{n}\equiv 0$ for all $n$ in (10) and this field remains at rest.
Looking at the equation for $\dot{p}_{n}$, either in its generic form (10) or
in its reduction to the D.N. shell model (11), we see that all but the viscous
term depend linearly on shell velocities $u_{n}$ and will therefore cancel out
if $u$ remains at rest. In this way, the Chernykh-Stepanov will converge to a
trivial instanton solution.

Therefore, one must take care to choose at least global forces or observables,
if not both. Moreover, local forces present ill-conditioned inverse
correlators. To sidestep both issues, we studied in detail how the system
responds to power law forcing correlations

$\chi_{nm}=k_{n}^{x}\delta_{nm}$

and choose a value of the exponent $x$ small enough ($x\leq-6$) such as to
mimic an almost local forcing as in realistic three dimensional turbulent
flows. On the opposite, for large exponent values, the system is fully driven
by the Gaussian forcing at all scales; this constitutes another interesting
limiting scenario for validation.

###  4.2 Instanton solutions

In this first work we focus on results about the simplest linear global
observable

$O[u]=\sum_{n}u_{n}.$ (12)

For each multiplier value $\lambda$, we seek a numerical approximation to the
corresponding instanton using the Chernykh-Stepanov iterative method. To this
end, we start with null fields $u$ and $p$ on a uniform grid on the time
interval $[-T,0]$. Then, we perform the iterative method described in
Subsection 3.3 while employing a suitable integrator; here, we used the two-
step Adam-Bashforth method. Figure 5 depicts an instanton solution obtained in
this way for $\lambda=2$, having converged after 175 interations. We observe
that as $t$ increases, more shells are excited and contribute to the
evaluation of the observable for which the instanton was defined; this can be
seen in the evolving time slices of the energy spectrum shown in Figure 6.

Figure 5: Instanton solution for (11) for $\lambda=2$,
$O=\sum_{n}u_{n}=82.34$; $c_{1}=10^{-3}$, $\nu=10^{-2}$, $c_{2}=0$ and
$\chi_{nm}=k_{n}^{-6}\delta_{nm}$. Figure 6: Squared velocities of each shell
in semi-log scale for an instanton approximation as it develops over different
time slices. Parameters used: $\lambda=2$, $O=\sum_{n}u_{n}=82.34$;
$c_{1}=10^{-3}$, $\nu=10^{-2}$, $c_{2}=0$ and
$\chi_{nm}=k_{n}^{-6}\delta_{nm}$. Figure 7: Squared velocities of different
shells for the final fully developed instantons ($t=0$) obtained for various
values of $\lambda$; curves are labelled according to the corresponding value
of the observable $O=\sum_{n}u_{n}$. Parameters used: $c_{1}=10^{-3}$,
$c_{2}=0$, $\nu=10^{-2}$ and $\chi_{nm}=k_{n}^{-6}\delta_{nm}$.

Having solved the instanton equations for a given value of $\lambda$, we
evaluate the observable and the action mapped to this value. By collecting
these values for a sufficient range of $\lambda$’s, we can evaluate the PDF
and the mean value of the observable using equations (6).

The validation of the methodology and implementation code was performed by
observing that in the fully linear case ($c_{1}=c_{2}=0$) each equation of the
shell becomes uncoupled and becomes an Ornstein-Uhlenbeck process. In this
setting, we can compare the distributions of an observable obtained from
direct numerical simulation (DNS) and instanton calculations, we also have
analytical results for the average and variance of the distributions of each
shell velocity.

###  4.3 Tail distributions for extreme fluctuations in the shell model

As one can see from Figure 7, instantons leading to high values of $O$ are
characterized by a non trivial energy distribution among many shells,
supporting the idea that the system is developing non trivial multi-scale
correlation to match the constrain. The PDF for an observable is computed up
to a prefactor from evaluation of the constrained action $S_{\lambda}$ for a
range of values of $\lambda$ as $\exp(-S_{\lambda})$. In order to understand
how well this PDF constructed only from instantons approximates the
distribution of the observable $O$ originally obtained from direct numerical
simulation of the shell model, we systematically explore the domain of
application of the methodology, starting with the fully linear setting we used
for validation and gradually increasing the values of the coefficient $c_{1}$.

In Figures 8 and 9 we show the PDFs of the observable compared with the one
from the numerical simulation (DNS) both in lin-lin and log-lin scale for two
values of non-linearity $c_{1}=10^{-4}$ and $10^{-3}$, respectively. As one
can see, the agreement is pretty good, not only for very intense values of the
observable. By increasing the non-linearity, the agreement deteriorate, as
shown in Figure 10 by the calculation of the excess kurtosis:

$F=\frac{\langle(O-\langle O\rangle)^{4}\rangle}{\langle(O-\langle
O\rangle)^{2}\rangle^{2}}-3.$

despite this, the tail of the distributions still shows good agreement also
for this relatively larger values of nonlinearities, as shown in Figure 9. By
increasing even further $c_{1}$, the agreement between the calculations based
only on the contributions from instantons and the full DNS is worsening.

Figure 8: Comparison of PDFs for the observable $O=\sum_{n}u_{n}$ obtained
from DNS and instanton calculations, linear scale to the left and log scale in
the y-axis to the right. Parameters used: $c_{1}=10^{-4}$, $c_{2}=0$,
$\nu=10^{-2}$ and $\chi_{nm}=k_{n}^{-6}\delta_{nm}$. Figure 9: Comparison of
PDFs for the observable $O=\sum_{n}u_{n}$ obtained from DNS and instanton
calculations, linear scale to the left and log scale in the y-axis to the
right. Parameters used: $c_{1}=10^{-3}$, $c_{2}=0$, $\nu=10^{-2}$ and
$\chi_{nm}=k_{n}^{-6}\delta_{nm}$. Figure 10: Comparison of excess kurtosis
for PDFs of the observable $O=\sum_{n}u_{n}$ for varying values of $c_{1}$
with $\nu=10^{-2}$, $c_{2}=0$ and $\chi_{nm}=k_{n}^{-6}\delta_{nm}$.

##  5 Discussion and conclusions

We conclude that the shell models, and especially the Desniansky-Novikov shell
model, are excellent candidates for testing the importance sampling
methodology, allowing for extensive sampling even of rare events with relative
simplicity if compared for full Navier-Stokes equations or other fluid models.
Starting from the ’Heat equations’ limit provides excellent validation ground
for the new methods, as the models reduce to a sequence of Ornstein-Uhlenbeck
processes, where analytical and numerical results are readily available. The
MSRJD formalism was applied to Shell Models in general, and instanton
equations were successfully defined and solved using the Chernykh-Stepanov
iterative method. We observed a good agreement between DNS and instanton-based
PDFs in the regime of weak nonlinearity where dynamics in all shells is still
mainly dominated by dissipative terms, even as the PDFs deviate from
Gaussianity due to ever increasing asymmetry in the distributions. However, as
soon as nonlinearity becomes relevant in the first few shells and outliers
increase in magnitude in the energy dissipation rate, we observe a difference
between DNS and instanton-based results. Instantons predict larger kurtosis
values and heavier positive tails. We expect that the results can be improved
by enhancing the method to take into account second-order fluctuations around
the instantons [41, 42, 43], here obtained by considering only saddle-point
approximations in the functional series expansion of the action. Improvements
in similar settings have been found by [11] for Burgers turbulence. Work in
this direction will be reported elsewhere. A very important question is also
connected to the set of observable used to constraint the action, and more
work must be done to explore instantons for small-scales quantities as the
gradients $\sum_{n}k_{n}u_{n}$ or the energy dissipation
$\sum_{n}k_{n}^{2}u_{n}^{2}$.

Extensions to D.N. models with complex variables, which presents more
interesting dynamics with chaotic attractors, would also be of interest. We
note that the formalism presented here is defined for shell models in general
by simply changing the functional $G_{n}[u]$ in the action (9) and instanton
equations (10). Instanton-based importance sampling in the GOY and SABRA shell
models for hydrodynamic turbulence [44, 45, 24] as well as in shell models
modeling magnetohydrodynamics [15, 16, 17], passive scalars [18, 19, 20] and
convective turbulence [21, 22] can also be investigated.

Acknowledgments We warmly thank Prof. Rainer Grauer and Timo Schorlepp (Ruhr-
University Bochum) for the valuable discussions and advice. We also thank
Prof. Karl Jansen (DESY-Zeuthen) for his valuable advice and for access to the
Warp cluster where simulations were performed. This work has received funding
from the European Union Horizon 2020 research and innovation programme under
the Marie Skłodowska-Curie grant agreement No 765048 and by the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (grant agreement No. 882340).

## References

  *   * Frisch [1995] Frisch, U.: Turbulence: the Legacy of AN Kolmogorov. Cambridge University Press, New York (1995)
  * Ohkitani and Yamada [1989] Ohkitani, K., Yamada, M.: Temporal intermittency in the energy cascade process and local lyapunov analysis in fully-developed model turbulence. Progress of theoretical physics 81(2), 329–341 (1989) https://doi.org/10.1143/PTP.81.329
  * Bohr et al. [1998] Bohr, T., Jensen, M.H., Paladin, G., Vulpiani, A.: Dynamical Systems Approach to Turbulence. Cambridge Nonlinear Science Series. Cambridge University Press, New York (1998). https://doi.org/10.1017/CBO9780511599972
  * Biferale [2003] Biferale, L.: Shell Models of Energy Cascade in Turbulence. Annual Review of Fluid Mechanics 35(1), 441–468 (2003) https://doi.org/10.1146/annurev.fluid.35.101101.161122
  * Dombre and Gilson [1998] Dombre, T., Gilson, J.L.: Intermittency, chaos and singular fluctuations in the mixed Obukhov-Novikov shell model of turbulence. Physica D: Nonlinear Phenomena 111(1-4), 265–287 (1998) https://doi.org/10.1016/S0167-2789(97)80015-2 arXiv:9510009 [chao-dyn]
  * Mahrt [1989] Mahrt, L.: Intermittency of atmospheric turbulence. J. Atmos. Sci. 46(1), 79–95 (1989) https://doi.org/10.1175/1520-0469(1989)046<0079:IOAT>2.0.CO;2"
  * Klymak et al. [2008] Klymak, J.M., Pinkel, R., Rainville, L.: Direct breaking of the internal tide near topography: Kaena ridge, hawaii. J. Phys. Oceanogr. 38(2), 380–399 (2008) https://doi.org/10.1175/2007JPO3728.1
  * Falkovich and Pumir [2008] Falkovich, G., Pumir, A.: Sling effect in collisions of water droplets in turbulent clouds. J. Atmos. Sci. 64(12), 4497–4505 (2008) https://doi.org/10.1175/2007JAS2371.1
  * Grafke et al. [2015] Grafke, T., Grauer, R., Schäfer, T., Vanden-Eijnden, E.: Relevance of instantons in Burgers turbulence. EPL (Europhysics Letters) 109(3), 34003 (2015) https://doi.org/10.1209/0295-5075/109/34003
  * Margazoglou et al. [2019] Margazoglou, G., Biferale, L., Grauer, R., Jansen, K., Mesterházy, D., Rosenow, T., Tripiccione, R.: Hybrid Monte Carlo algorithm for sampling rare events in space-time histories of stochastic fields. Physical Review E 99(5), 053303 (2019) https://doi.org/10.1103/PhysRevE.99.053303 arXiv:1808.02020
  * Apolinário et al. [2019] Apolinário, G., Moriconi, L., Pereira, R.: Onset of intermittency in stochastic burgers hydrodynamics. Physical Review E 99(3), 033104 (2019) https://doi.org/10.1103/PhysRevE.99.033104
  * Grafke et al. [2013] Grafke, T., Grauer, R., Schäfer, T.: Instanton filtering for the stochastic burgers equation. Journal of Physics A: Mathematical and Theoretical 46(6), 062002 (2013) https://doi.org/10.1088/1751-8113/46/6/062002
  * Migdal [2020] Migdal, A.: Clebsch confinement and instantons in turbulence. International Journal of Modern Physics A 35(31), 2030018 (2020) https://doi.org/10.1142/S0217751X20300185
  * Moriconi and Pereira [2022] Moriconi, L., Pereira, R.M.: Statistics of extreme turbulent circulation events from multifractality breaking. Phys. Rev. E 106, 054121 (2022) https://doi.org/%****␣main.tex␣Line␣675␣****10.1103/PhysRevE.106.054121
  * Frick and Sokoloff [1998] Frick, P., Sokoloff, D.: Cascade and dynamo action in a shell model of magnetohydrodynamic turbulence. Physical Review E 57(4), 4155 (1998) https://doi.org/10.1103/PhysRevE.57.4155
  * Nigro et al. [2004] Nigro, G., Malara, F., Carbone, V., Veltri, P.: Nanoflares and mhd turbulence in coronal loops: a hybrid shell model. Physical Review Letters 92(19), 194501 (2004) https://doi.org/10.1103/PhysRevLett.92.194501
  * Plunian et al. [2013] Plunian, F., Stepanov, R., Frick, P.: Shell models of magnetohydrodynamic turbulence. Physics Reports 523(1), 1–60 (2013) https://doi.org/10.1016/j.physrep.2012.09.001
  * Jensen et al. [1992] Jensen, M., Paladin, G., Vulpiani, A.: Shell model for turbulent advection of passive-scalar fields. Physical Review A 45(10), 7214 (1992) https://doi.org/10.1103/PhysRevA.45.7214
  * Cohen et al. [2002] Cohen, Y., Gilbert, T., Procaccia, I.: Statistically preserved structures in shell models of passive scalar advection. Physical Review E 65(2), 026314 (2002) https://doi.org/10.1103/PhysRevE.65.026314
  * Mitra and Pandit [2005] Mitra, D., Pandit, R.: Dynamics of passive-scalar turbulence. Physical review letters 95(14), 144501 (2005) https://doi.org/10.1103/PhysRevLett.95.144501
  * Brandenburg [1992] Brandenburg, A.: Energy spectra in a model for convective turbulence. Physical review letters 69(4), 605 (1992) https://doi.org/10.1103/PhysRevLett.69.605
  * Ching [2010] Ching, E.S.: Studying anomalous scaling and heat transport of turbulent thermal convection using a dynamical model. Physica D: Nonlinear Phenomena 239(14), 1346–1352 (2010) https://doi.org/10.1016/j.physd.2009.10.021
  * Campolina and Mailybaev [2018] Campolina, C.S., Mailybaev, A.A.: Chaotic blowup in the 3d incompressible euler equations on a logarithmic lattice. Physical review letters 121(6), 064501 (2018)
  * Ditlevsen [2011] Ditlevsen, P.: Turbulence and Shell Models. Cambridge University Press, New York (2011). https://doi.org/10.1017/CBO9780511919251
  * Düben et al. [2008] Düben, P., Homeier, D., Jansen, K., Mesterhazy, D., Münster, G., Urbach, C.: Monte Carlo simulations of the randomly forced Burgers equation. EPL (Europhysics Letters) 84(4), 40002 (2008) https://doi.org/10.1209/0295-5075/84/40002
  * Farazmand and Sapsis [2017] Farazmand, M., Sapsis, T.P.: A variational approach to probing extreme events in turbulent dynamical systems. Science Advances 3(9), 1701533 (2017) https://doi.org/10.1126/sciadv.1701533
  * Mesterházy and Jansen [2011] Mesterházy, D., Jansen, K.: Anomalous scaling in the random-force-driven Burgers’ equation: a Monte Carlo study. New Journal of Physics 13(10), 103028 (2011) https://doi.org/10.1088/1367-2630/13/10/103028
  * Schafer [2004] Schafer, T.: Instantons and Monte Carlo methods in quantum mechanics (2004) https://doi.org/10.48550/arXiv.hep-lat/0411010
  * Martin et al. [1973] Martin, P.C., Siggia, E.D., Rose, H.A.: Statistical Dynamics of Classical Systems. Physical Review A 8(1), 423–437 (1973) https://doi.org/10.1103/PhysRevA.8.423
  * Janssen [1976] Janssen, H.-K.: On a Lagrangean for classical field dynamics and renormalization group calculations of dynamical critical properties. Zeitschrift für Physik B Condensed Matter and Quanta 23(4), 377–380 (1976) https://doi.org/10.1007/BF01316547
  * Dominicis [1976] Dominicis, C.d.: Technics of field renormalization and dynamics of critical phenomena. In: J. Phys.(Paris), Colloq, pp. 1–247 (1976). https://doi.org/10.1051/jphyscol:1976138
  * Onsager and Machlup [1953] Onsager, L., Machlup, S.: Fluctuations and Irreversible Processes. Physical Review 91(6), 1505–1512 (1953) https://doi.org/10.1103/PhysRev.91.1505
  * Grafke et al. [2015] Grafke, T., Grauer, R., Schäfer, T.: The instanton method and its numerical implementation in fluid mechanics. Journal of Physics A: Mathematical and Theoretical 48(33), 333001 (2015) https://doi.org/%****␣main.tex␣Line␣950␣****10.1088/1751-8113/48/33/333001
  * Stratonovich [1957] Stratonovich, R.: On a method of calculating quantum distribution functions. In: Soviet Physics Doklady, vol. 2, p. 416 (1957)
  * Hubbard [1959] Hubbard, J.: Calculation of partition functions. Physical Review Letters 3(2), 77 (1959)
  * Duane et al. [1987] Duane, S., Kennedy, A.D., Pendleton, B.J., Roweth, D.: Hybrid monte carlo. Physics Letters B 195(2), 216–222 (1987) https://doi.org/10.1016/0370-2693(87)91197-X
  * Lüscher [2010] Lüscher, M.: Computational Strategies in Lattice QCD (2010). https://doi.org/%****␣main.tex␣Line␣1000␣****10.48550/arXiv.1002.4232
  * Brooks et al. [2011] Brooks, S., Gelman, A., Jones, G., Meng, X.-L. (eds.): Handbook of Markov Chain Monte Carlo. Chapman and Hall/CRC, Florida (2011). https://doi.org/10.1201/b10905
  * Chernykh and Stepanov [2001] Chernykh, A.I., Stepanov, M.G.: Large negative velocity gradients in burgers turbulence. Phys. Rev. E 64, 026306 (2001) https://doi.org/10.1103/PhysRevE.64.026306
  * Ebener et al. [2019] Ebener, L., Margazoglou, G., Friedrich, J., Biferale, L., Grauer, R.: Instanton based importance sampling for rare events in stochastic PDEs. Chaos 29 (2019) https://doi.org/10.1063/1.5085119
  * Schorlepp et al. [2021] Schorlepp, T., Grafke, T., Grauer, R.: Gel’fand–yaglom type equations for calculating fluctuations around instantons in stochastic systems. Journal of Physics A: Mathematical and Theoretical 54(23), 235003 (2021) https://doi.org/10.1088/1751-8121/abfb26
  * Moriconi et al. [2014] Moriconi, L., Pereira, R., Grigorio, L.: Velocity-gradient probability distribution functions in a lagrangian model of turbulence. Journal of Statistical Mechanics: Theory and Experiment 2014(10), 10015 (2014) https://doi.org/10.1088/1742-5468/2014/10/P10015
  * Apolinário et al. [2019] Apolinário, G., Moriconi, L., Pereira, R.: Instantons and fluctuations in a lagrangian model of turbulence. Physica A: Statistical Mechanics and its Applications 514, 741–757 (2019) https://doi.org/10.1016/j.physa.2018.09.102
  * Yamada and Ohkitani [1988] Yamada, M., Ohkitani, K.: Lyapunov spectrum of a model of two-dimensional turbulence. Phys. Rev. Lett. 60, 983–986 (1988) https://doi.org/10.1103/PhysRevLett.60.983
  * L’vov et al. [1998] L’vov, V.S., Podivilov, E., Pomyalov, A., Procaccia, I., Vandembroucq, D.: Improved shell model of turbulence. Phys. Rev. E 58, 1811–1822 (1998) https://doi.org/10.1103/PhysRevE.58.1811


# Impact of tropical teleconnections on the long-range predictability of the
atmosphere at midlatitudes: A reduced-order multi-scale model perspective

Stéphane Vannitsem  Royal Meteorological Institute of Belgium, Avenue
Circulaire, 3, 1180 Brussels, Belgium

###### Abstract

Teleconnections between the tropical and the extratropical climates are often
considered as a potential source of long-term predictability at seasonal to
decadal time scales in the extratropics. This claim is taken up in the present
work by investigating the predictability of a coupled ocean-atmosphere
extratropical model under a one-way forcing generated by a tropical model.
Both models display a chaotic dynamics, and the dominant variable of the
extratropical model displays a high correlation with the tropical forcing in
the reference simulation, inducing a low-frequency variability signal in the
extratropics.

Numerical experiments emulating the presence of initial condition errors in
the tropical model are conducted to clarify their impact on the predictability
in the extratropics. It is shown that: (i) The correlation skill of the
dominant observable affected by the forcing is considerably degraded at
interannual time scales due to the presence of initial condition errors in the
tropics, considerably limiting the potential of teleconnections; (ii)
averaging of an ensemble of forecasts – with a small number of members – may
substantially improve the quality of the forecasts; and (iii) temporal
averaging may also improve the quality of the forecasts (at the expense of
being able to forecast extreme events), in particular when the forcing affects
weakly the observable under interest.

<EMAIL_ADDRESS>

  * •

May 20, 2023

Keywords: Teleconnections, low-frequency variability, ENSO, predictability,
chaos

##  1 Introduction

The tropical regions display key low-frequency variability ranging from months
with the Madden-Julian Oscillation to annual and decadal time scales for the
El-Niño-Southern Oscillation (ENSO). These processes display teleconnections
with large regions in the extratropics as discussed in [1, 2, 3, 4, 5]. These
are often expressed in terms of correlations between the tropical and
extratropical regions, suggesting a potential influence on one another. Note
however than correlation does not mean causation as already stressed in many
recent works on causality analyses [6, 7, 8, 9]. In [9], a new approach of
causality detection was used to clarify the transfer of information among
different climate indices and in particular ENSO on monthly to decadal time
scales. It was in particular shown that ENSO has considerable influence on the
Pacific north Atlantic pattern, on the Pacific decadal oscillation, on the
tropical north Atlantic region and slightly on the Arctic oscillation. The
mechanism often put forward for the teleconnections with the midlatitudes is
the atmospheric bridge described in [2].

As these tropical processes have typical time scales of evolution much longer
than the ones of the extratropical dynamics, it is natural to expect that
these low-frequency variability signals would emerge in the evolution of the
extratropical atmosphere, and could be exploited in our seasonal and decadal
forecasts of the midlatitude dynamics [10, 11], although the way these
frequencies will emerge is a complicate matter as illustrated in [12]. But the
processes developing in the tropical regions are not simple periodic signals
that could be easily exploited as for instance the annual cycle. They display
an erratic behavior reminiscent of chaotic or stochastically-forced dynamics
[13]. This implies that these processes display the property of sensitivity to
initial conditions that will affect any forecast of the tropical evolution as
initial state errors are always present in the model forecasts, and these
uncertainties will, in turn, affect the evolution of the signal at
midlatitudes.

The starting point of the present work is the claim that any forecast will
hardly be able to exploit the full potential of teleconnections. This claim is
demonstrated by investigating the impact of errors in the initial conditions
of a tropical forcing on an extratropical model. This system has been recently
built to investigate the emergence of pullback attractors in a midlatitude
reduced-order coupled ocean-atmosphere model forced by a tropical low-order
model [12]. The midlatitude model is the so-called VDDG model (for Vannitsem-
Demaeyer-DeCruz-Ghil) describing the interaction between a 2-layer quasi-
geostrophic atmosphere with a 1.5-layer quasi-geostrophic ocean [14], while
the tropical low-order model is a box model developed by [15, 16, 13]. A one-
way interaction is imposed from the tropical model to the extratropical one.
The solutions of both models that will be used in the present work display
chaotic dynamics, and are therefore sensitive to the initial conditions. The
experiments that will be performed will be to introduce random initial state
errors of different amplitudes in the tropical model to clarify the ability of
the extratropical model to produce long-term skillful forecasts at
midlatitudes.

In Section 2, the coupled tropical-extratropical model setup is described,
together with the parameter used in the current experiment and the reference
solutions around which the predictability experiments will be performed.
Section 3 is devoted to the description of the experimental setup and the
analysis of the different experiments. The impact of ensemble averaging and
temporal averaging are also explored in this section. The main conclusions are
drawn in Section 4.

##  2 Governing equations for the tropical–extratropical model

###  2.1 The ENSO module

The ENSO model was developed in a series of papers by F.-F. Jin, A. Timmermann
and colleagues [15, 16, 17, 13, 18]. They modeled the dynamics of the ocean’s
upper layer in the tropical Pacific using a low number of variables. Their
ENSO model describes the dynamics of the temperature in the eastern and
western tropical Pacific basins, completed by an equation for the evolution of
the thermocline depth. The model represents the horizontal discharge-recharge
mechanisms at play in the tropical Pacific through the heat exchanges between
the tropical and subtropical waters, subject to surface wind stress and
upwelling of subsurface cold water in the eastern part of the domain [16].
[18] developed a nondimensional model version in which the time, for instance,
is normalized by a typical time scale of tropical wave propagation of roughly
3.5 months. The latter ENSO model version is coupled here with an
extratropical module whose time is nondimentionalized by the Coriolis
parameter $f_{0}$, and the [18] equations are therefore slightly modified in
order to incorporate this change of reference time scale:

$\displaystyle\frac{dx}{dt}$
$\displaystyle=\rho\delta(x^{2}-ax)+sx(x+y+c-c\tanh(x+z)),$ (1a)
$\displaystyle\frac{dy}{dt}$ $\displaystyle=-\rho\delta(ay+x^{2}),$ (1b)
$\displaystyle\frac{dz}{dt}$ $\displaystyle=\delta(k-z-\frac{x}{2}).$ (1c)

Here $x,y,z$ and $t$ are dimensionless, $x$ is the temperature difference
between the eastern and western basins of the tropical Pacific, $y$ the
western basin’s temperature anomaly with respect to a reference value, and $z$
the western basin’s thermocline depth anomaly.

The dimensionless parameters are defined as follows

$\displaystyle a$ $\displaystyle=\frac{\alpha bL}{\epsilon
h^{\\*}\beta},\qquad\rho=\frac{\epsilon
h^{\\*}\beta}{rbL},\qquad\delta=\frac{r}{f_{0}},$ (2a) $\displaystyle c$
$\displaystyle=\frac{C}{S_{0}},\qquad k=\frac{K}{S_{0}},\qquad s=\frac{\zeta
h^{\\*}\beta}{bLf_{0}},$ (2b)

with $f_{0}$ the Coriolis parameter used to nondimensionalize the time, and
the other parameters as defined in [18, Table 1]. For the applications
presented in this paper, one set of parameter values is considered leading to
a chaotic solution, see Table 1.

Table 1: Dimensionless parameter values of the ENSO module $a=7.658609809$
---
$\rho=0.29016$
$\delta=0.0002803$
$c=2.3952$
$k=0.4032$
$s=0.001069075$

The ENSO module is taken here to be unaffected by the midlatitude module, and
is thus the driving system in the coupled tropical-extratropical model.
Equations (1) can thus be integrated independently of the rest of the model,
and this is done using a fourth-order Runge-Kutta scheme with a time step
$\Delta t=0.1346$ hours $=0.05$ nondimensional time units. The variable $x+y$
corresponds in the model (1) to the sea surface temperatures in the eastern
tropical Pacific that are commonly associated with the Niño-3 index, that will
be referred later as $T_{ENSO}$. A long chaotic solution was started from the
initial state $(x=-2.8439,y=-0.62,z=1.480)$ that will be used in the following
as the long reference trajectory of the tropical dynamics.

###  2.2 The VDDG extratropical model

The coupled ocean–atmosphere model used herein for the midlatitudes was
developed by [19] and it is called hereafter the VDDG model. Different
versions of this model have been used to study low-frequency variabiliy (LFV)
within the coupled ocean–atmosphere system [20, 21], to investigate the
properties of the Lyapunov exponents in such a system [22, 23], to build the
stochastic parametrization of subgrid-scale forcing [24, 25], and to develop
data assimilation schemes in coupled models [26, 27, 28].

The atmospheric module is based on the vorticity equations of a two-layer
quasi-geostrophic flow defined on a beta-plane [29, 30]. The equations in
pressure coordinates are:

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

Here $\psi^{1}_{a}$ and $\psi^{3}_{a}$ are the streamfunction fields at 250
hPa and 750 hPa, respectively, while $\omega=dp/dt$ is the vertical velocity,
$f_{0}$ the Coriolis parameter at $\phi_{0}=$ 45 ∘ latitude, and $\beta=df/dy$
at $\phi_{0}$. The coefficients $k_{d}$ and $k^{\prime}_{d}$ multiply the
surface friction term and the internal friction between the layers,
respectively.

An additional term has been introduced in Eq. (3b) in order to account for the
presence of a surface boundary velocity $\psi_{o}$ of the oceanic flow; see
Eq. (4) below. This term corresponds to the Ekman pumping on a moving surface
and is the mechanical contribution of the interaction between the ocean and
the atmosphere.

The ocean module is described by the reduced-gravity, quasi-geostrophic
shallow-water model [29, 30]. The forcing is provided by the wind generated by
the atmospheric module above. The governing equation is:

$\frac{\partial}{\partial
t}\left(\nabla^{2}\psi_{o}-\frac{\psi_{o}}{L_{R}^{2}}\right)+J(\psi_{o},\nabla^{2}\psi_{o})+\beta\frac{\partial\psi_{o}}{\partial
x}=-r\nabla^{2}\psi_{o}+\frac{{\mathrm{curl}}_{z}\vec{\tau}}{\rho H},$ (4)

where $\psi_{o}$ is the streamfunction, $\rho$ the density of water, $H$ the
depth of the fluid layer, $L_{R}$ the reduced Rossby deformation radius, $r$ a
Rayleigh friction coefficient at the bottom of the fluid layer, and
curl${}_{z}\vec{\tau}$ is the vertical component of the wind stress curl.

The wind stress in the VDDG model is given by
$(\tau_{x},\tau_{y})=C(u-U,v-V)$, where $(u=-\partial\psi^{3}_{\rm a}/\partial
y,v=\partial\psi^{3}_{\rm a}/\partial x)$ are the horizontal components of the
geostrophic wind, and $(U,V)$ the components of the geostrophic currents in
the ocean. One thus gets

$\mathrm{curl}_{z}\tau=C\nabla^{2}(\psi^{3}_{\rm a}-\psi_{\rm o}).$ (5)

and the wind stress is proportional to the relative velocity between the flow
in the ocean’s upper layer and the wind in the lower atmospheric layer. The
drag coefficient $d=C/(\rho_{0}H)$ gives the strength of the mechanical
coupling between the ocean and the atmosphere and it was a key bifurcation
parameter in [19].

The dynamic equations (3) and (4) are supplemented by temperature equations
for the two subsystems. For the atmosphere,

$\gamma_{a}\left(\frac{\partial T_{a}}{\partial
t}+J(\psi_{a},T_{a})-\sigma\omega\frac{p}{R}\right)=-\lambda(T_{a}-T_{o})+E_{a,R},$
(6)

with

$E_{a,R}=\epsilon_{a}\sigma_{B}T_{o}^{4}-2\epsilon_{a}\sigma_{B}T_{a}^{4}+R_{a}.$
(7)

Here $R$ is the gas constant, $\epsilon_{a}$ the emissivity of the atmosphere,
$\sigma_{B}$ the Stefan-Boltzman constant, $R_{a}$ the shortwave radiation at
the top of the atmosphere, $\omega$ the vertical velocity in pressure
coordinates, and $\sigma=-(R/p)(\partial T_{a}/\partial p-1/(\rho_{a}c_{p}))$
is the static stability, with $p$ the pressure, $\rho_{a}$ the air density,
$c_{p}$ the specific heat at constant pressure, and $\sigma$ here is taken to
be a constant. Note also that, thanks to the hydrostatic relation in pressure
coordinates and the ideal gas relation $p=\rho_{a}RT_{a}$, the atmospheric
temperature can be written as $T_{a}=-(p/R)f_{0}(\partial\psi_{a}/\partial
p)$.

Table 2: List of parameters of the extratropical VDDG module Parameter (unit) | Value | Parameter (unit) | Value
---|---|---|---
$L_{y}=\pi L$ (km) | $5.0\times 10^{3}$ |  $\gamma_{\rm{o}}$ (J m-2 K-1) | $4\times 10^{6}\,h$
$f_{0}$ (s-1) | $1.032\times 10^{-4}$ |  $C_{\rm{o}}$ (W m-2) | 310
$n=2L_{y}/L_{x}$ | $1.5$ |  $T_{\rm{o}}^{0}$ (K) | $285$
$R_{\rm{E}}$ (km) | $6370$ |  $\gamma_{\rm{a}}$ (J m-2 K-1) | $1.0\times 10^{7}$
$\phi_{0}$ | $\pi/4$ |  $C_{\rm{a}}$ (W m-2) | $C_{\rm{o}}/4$
$g^{\prime}$ | $3.1\times 10^{-2}$ | $\epsilon_{\rm{a}}$ | $0.76$
$r$ (s-1) | $1.0\times 10^{-7}$ |  $\beta$ (m-1 s-1) | $1.62\times 10^{-11}$
$H$ (m) | $100$ |  $T_{\rm{a}}^{0}$ (K) | $270$
$d$ (s-1) | $C/(\rho_{o}H)$ |  $\lambda$ (W m-2 K-1) | $1004\,C$
$k_{d}$ (s-1) | $(gC)/(\Delta p)$ |  $R$ (J kg-1 K-1) | $287$
$k_{d}^{\prime}$ (s-1) | $(gC)/(\Delta p)$ |  $\sigma$ (J kg-1 Pa-2) | $2.16\times 10^{-6}$
$C$ (kg m-2 s-1) | $0.008$ |  |

For the ocean,

$\gamma_{o}\left(\frac{\partial T_{o}}{\partial
t}+J(\psi_{o},T_{o})\right)=-\lambda(T_{o}-T_{a})+E_{R},$ (8)

with

$E_{R}=-\sigma_{B}T_{o}^{4}+\epsilon_{a}\sigma_{B}T_{a}^{4}+R_{o}.$ (9)

Here $R_{o}$ is the shortwave radiation entering the ocean, $\gamma_{o}$ the
heat capacity of the ocean, and $\lambda$ is the inverse of the time scale
associated with heat transfer between the ocean and the atmosphere, which
includes both the latent and sensible heat fluxes. In fact, we assume that
this combined heat transfer is proportional to the temperature difference
between the atmosphere and the ocean.

The temperatures in both modules are linearized around a reference value in
order to reduce the quartic terms of the energy balance equations (7) and (9)
to linear terms, assuming that the temperature fluctuations are small. This
modification helps one to reduce the number of terms on the right-hand side of
the ordinary differential equations (6) and (8) when building the spectral
low-order model.

The model fields in both its atmosphere and its ocean are developed in Fourier
series and truncated at a low order. The number of modes herein is fixed at 10
for the atmosphere and 8 for the ocean, leading to 20 ordinary differential
equations for the former and 16 for the latter. This model configuration is
the original VDDG one; see also [21]. The parameter values used in the present
work are listed in Table 2. In the configuration of Table 2, there is no
substantial low-frequency variability emerging in the extratropical component
(See [12]), allowing for solely clarifying the impact of external low-
frequency variability.

###  2.3 Modeling the tropical–extratropical interaction

[31] investigated the changes of extratropical wintertime cyclogenesis when El
Niño or La Niña events are occurring in the tropical Pacific. They showed, in
particular, that the background zonal-flow anomaly is more intense over the
north Atlantic during La Niña, while it is stronger over the north Pacific
during El Niño. This finding tells us that an important effect of the tropical
forcing is to change the intensity of the zonal flow in either region, and
that the impact of El Niño and La Niña differs from one region to the other.
In order to mimic this dynamic effect in the extratropical VDDG model used
here, we impose a direct linear forcing of the model’s first barotropic
atmospheric mode. It is this barotropic streamfunction mode that represents
the intensity of the zonal flow within the atmosphere. Its dynamics is written
as

$\frac{d\psi_{a,1}}{dt}=f_{1}(\psi_{a,1},\theta_{a,1})+g\delta(x+y).$ (10)

Here $f_{1}(\psi_{a,1},\theta_{a,1})$ is the original right-hand side of the
dynamical evolution (3a) of $\psi_{a,1}$; $\delta(x+y)$ represents the eastern
tropical Pacific basin’s temperature anomalies; and $g$ scales the intensity
of the tropical forcing.

Thus $g$ represents the crucial forcing of the midlatitude VDDG model
described in Sec. 2.2 above by the ENSO module of Sec. 2.1. In our setting,
given a positive $g$-value, a positive, warm anomaly will induce an increase
of $\psi_{a,1}$, and hence of the mean zonal flow,
$U=-\psi_{a,1}\partial(\sqrt{2}\cos y)/\partial y$. This situation corresponds
to the intensification of the zonal flow over the north Pacific during an El
Niño. If, to the contrary, $g$ is negative, this would correspond to an
intensification during La Niña that mimics the ENSO effect over the north
Atlantic. In the current analysis, we are focusing on one specific positive
value of $g$ providing a strong teleconnection between the tropics and the
midlatitudes.

Figure 1: Solutions of the coupled tropical-extratropical coupled ocean-
atmosphere model for the coupling parameter $g=0.002$. (a) The tropical
temperature forcing, $T_{ENSO}$; (b) the forced midlatitude variable,
$\psi_{a,1}$; (c) the first baroclinic mode, $\theta_{a,1}$; and (d) the
scatter plot of the variables $T_{ENSO}$ and $\psi_{a,1}$.

Figure 1 displays a typical trajectory of the idealized ENSO model for a
coupling parameter of $g=0.002$. Panel (a) shows the forcing trajectory of sea
surface temperature in the east equatorial pacific with a succession of strong
and less strong El-Niño, and la Niña events, occurring in an irregular manner.
The evolution of the first mode of the barotropic streamfunction,
$\psi_{a,1}$, is shown in panel (b), indicating regular bursting associated
with strong El-Niño events. Panel (c) shows the same type of evolution for the
first mode of the baroclinic streamfunction, $\theta_{a,1}$, strongly related
to the evolution of $\psi_{a,1}$. In this evolution, the bursting is not
visible anymore, indicating a much lower influence of the tropical forcing
than for $\psi_{a,1}$. Finally in panel (d), the relationship between the
instantaneous values of $\psi_{a,1}$ and the temperature in the east part of
the tropical ocean basin is displayed, indicating a strong teleconnection
between the modelled tropical Pacific and the streamflow of the midlatitude
model. The correlation between the two observables is equal to 0.914, a strong
teleconnection between the two regions in this idealized modelling. If the
tropical-extratropical coupling parameter $g$ is modified, one can either
increase or decrease the strength of these teleconnections. Such an impact
will be addressed in the next section.

In order to further clarify the link between the tropical and extratropical
systems, power spectra of $\psi_{a,1}$ and $T_{ENSO}$ for the coupling
parameter $g=0.002$ are displayed in Fig. 2. For the forcing $T_{ENSO}$, the
energy is concentrated at frequencies below 0.001, with a rapid drop at larger
frequencies. For $\psi_{a,1}$, the spectrum is much flatter at high
frequencies than for the forcing, indicating that a substantial amount of
energy is contained at that frequencies. Interestingly a strong similarity of
the spectra between frequencies of $10^{-5}$ and $0.002$ is found, indicating
the strong link between the forcing and $\psi_{a,1}$. Note that the broadband
peak of the $T_{ENSO}$ is obviously reflecting the chaotic nature of the
dynamics, but we can isolate a few dominating peaks at about 30, 15, 11 and 5
years. The latter is the one familiar in many time series analysis of the ENSO
evolution. The former ones may be emulating the long term behavior of ENSO
which displays stronger El-Niño events on decadal and multi-decadal time
scales [32, 33].

Figure 2: Power spectra of $\psi_{a,1}$ and $T_{ENSO}$ for the coupling
parameter $g=0.002$, as obtained using the Multi-Taper Method [34] with 50000
data points sampled every 2.804 days.

##  3 Predictability at midlatitudes

###  3.1 Experimental setup

The tropical-extratropical coupled system presented in section 2 allows for a
clean analysis of the limits of predictability of the extratropical climate
system forced by the tropics. As briefly mentioned in the previous section,
the solution generated by the tropical model is chaotic, implying that the
signal to be extracted also exhibits sensitivity to initial conditions as for
the extratropics. The Lyapunov time scale of the tropical solution is of the
order of 7 years, indicating a very slow increase of errors in time compatible
with the dynamics of ENSO.

A very long control run of the model is first generated which provides the
reference against which all forecasting experiments will be evaluated. This
constitutes our truth. For the forecasting experiments, it is assumed that an
estimate of the truth has been done at specific times along the reference
trajectory with some uncertainty. This uncertainty is generated using random
perturbations. We suppose that this uncertainty is coming from a uniform
distribution for the variables in both the tropical region and at
midlatitudes. It is also further assumed that we know this distribution
exactly and one can therefore generate an appropriate ensemble forecasting
system based on the same distribution as the one used to generate the
uncertainty of the observations. The ensemble forecasts are therefore done as
follows.

A first set of random perturbations mimicking the uncertainty on the initial
conditions in the extratropics is introduced. The distribution used is a
uniform distribution of mean 0 and range $[-5.10^{-7}:5.10^{-7}]$. This
amplitude is very small as compared with the amplitude of the variables (in
non-dimensional units). Typically of the order of 0.001 $\%$ of the
variability of $\psi_{a,1}$. This mimics a situation for which the initial
conditions are very well known in the extratropics, that will allow for having
a clear picture of the different phases of the error dynamics and the skill of
the model variables.

The tropical dynamics is also affected by initial condition errors, that will
in turn affect the forecasts in the extratropics. In order to mimic this
feature, a random initial error is also introduced in the initial conditions
of the tropical model, generated from a uniform distribution of zero mean and
an amplitude, $PT$, of either 0., 0.02 and 0.05, 0.10, 0.15 and 0.20. Thus the
distribution is covering the domain $[-PT/2:PT/2]$. The perturbation
amplitudes are therefore ranging from 0 to 15 $\%$ of the natural variability
of the tropical model temperature, $x+y$, used to force midlatitudes.

As in many long-term forecasting systems, an ensemble of forecasts is
generated starting from different initial conditions. The number, $M$, of such
ensemble members will be modified in order to clarify their impact on the
ensemble mean estimate. The initial conditions of these ensemble members are
perturbed consistently with the uncertainties affecting the two models in
order to generate a reliable forecast: Once an initial state is defined with
the uncertainties mentioned above, random perturbations sampled from the same
uniform distributions are generated and added to this initial state to produce
each member of the ensemble.

For verification, a few usual metrics are computed: the correlation between
the reference trajectory and the forecasts or the ensemble mean (referred as
the correlation skill), and the mean square error. To get estimates that are
independent of the initial conditions on the attractor of the coupled system,
N=2,000 realizations are performed starting at different location along the
reference trajectory.

###  3.2 Perfect Teleconnection Predictability

In order to illustrate the potential impact of the teleconnections introduced
in this system on the skill of the solutions, let us first consider the
experiments done without error in the initial conditions of the tropical
forcing for values of the coupling parameter $g=0$ and $0.002$. In other
words, it is assumed that the forcing evolution is perfectly known, which
should allow for exploiting optimally the teleconnections revealed in Fig. 1.
The error in the initial states of the extratropical model is fixed as
discussed in the previous section.

Figure 3a shows the correlation skill between the control integration (without
initial condition errors) and single forecasts ($M=1$) for variables
$\psi_{a,1}$ and $\theta_{a,1}$ for $g=0.$ and $g=0.002$. For $g=0.$, the
correlation skill is completely lost after about 50 days. When the forcing is
introduced with $g=0.002$, the picture is similar for the variable
$\theta_{a,1}$, while for variable $\psi_{a,1}$ the correlation skill is first
quickly decreasing and then saturates around a value of about 0.82. Two
general comments are in order here with the specific design of the tropical-
extratropical model: (i) The teleconnection present between the forcing and
$\psi_{a,1}$ induce a long term correlation skill of this variable with a
plateau at a level defined by the strength of the teleconnection induced by
the tropical forcing; and (ii) there is no obvious propagation of the
teleconnection dependence to the baroclinic variable, $\theta_{a,1}$, strongly
related to the dynamics of $\psi_{a,1}$. Although this type of consideration
is not entirely new, it provides a key reference for the experiments that will
be done next. In particular, the strong correlation between the reference and
forecasting trajectories for $\psi_{a,1}$ reveals the strong potential of
teleconnections. This correlation skill level for $g=0.002$ of $\psi_{a,1}$
will be referred in the following as the Perfect Teleconnection Predictability
(PTP).

Figure 3b shows the PTP for various values of $g$ from $0.$ to $0.003$. It
reveals that once $g$ is increased, the PTP increases accordingly. It is
however interesting to point to the fact that the level of correlation skill
shown from, say, $100$ to $500$ days, is always smaller than the correlation
between $\psi_{a,1}$ and $T_{ENSO}$: for $g=0.0005$, the correlation is
$0.50$, while the level of the corresponding correlation skill in Fig. 3b is
lower than $0.25$. This reflects the complex interplay between the forcing and
the forced variable. So the teleconnection as defined by the correlation
between $\psi_{a,1}$ and $T_{ENSO}$ is an upper bound of the actual
correlation skill of $\psi_{a,1}$.

From Fig. 3b, it is also interesting to note that whatever the value of $g$,
the quick initial decrease of skill toward the plateau is very similar. The
reason is that the short term predictability of the extratropical model is
almost the same when measured by the first Lyapunov exponent for this range of
values of $g$ (see [12]).

Figure 3: Correlation between the reference integration, viewed as perfect
observations, and the ensemble mean of forecasts for (a) variables $\psi_{1}$
and $\theta_{1}$ for $g=0.$ and $g=0.002$ and (b) variable $\psi_{1}$ for
$g=0.$, $0.001$, $0.002$ and $0.003$. The ensemble size is fixed to $M=1$,
hence corresponding to single deterministic forecasts.

###  3.3 Predictability for M=1

Let us now focus on the case $g=0.002$ and consider the more realistic
situation for which there are some uncertainties in the initial conditions of
the tropical forcing, but still with only one ensemble member $M=1$. As the
tropical model is chaotic, it also displays sensitivity to initial conditions
that will affect in turn the extratropical forecasts.

Figure 4 displays the evolution of the Mean Square Error (MSE) for several
amplitudes of the perturbations of the initial states of the tropical model.
On this doubly logarithmic scale, a first interesting information is the
quadratic increase of the MSE for short lead times. This behavior is well
documented and typical when model errors or boundary forcing errors are
present [35, 36, 37]. Interestingly this error quickly amplifies suggesting
that boundary forcing errors can considerably harm the forecasts in the
extratropics whatever the quality of the initial conditions of the
extratropical model. Another important information, is the level of saturation
of the error which is reached when forcing errors are present: The error at
long lead times is higher. Finally, the way it saturates from, say, 30 days to
the end of the forecasts shows a complex feature related to the time scales of
the forcing: A plateau is present between 30 to 100 days which is a period for
which the difference between the forcing acting on the reference and the
forecasting model is still small, and hence following the error evolution in
the absence of any initial condition error of the tropical model; and then
followed by an oscillating behavior associated with the typical oscillations
found in the original tropical model dynamics. The specific sharp decrease of
the error around 900 days needs further clarifications by investigating the
detailed evolution of the distribution of the errors as a function of time on
these time scales and beyond for various initial error amplitude in the
forcing. This question will be addressed in the future.

Figure 4: MSE evolution for the variable $\psi_{a,1}$ under the influence of
tropical forcing, affected by different initial condition error amplitudes,
$PT=0,0.02,0.05,0.10,0.15$.

Although the error evolution is an interesting quantity to analyze in
particular for the short terms, the correlation skill is probably more
relevant for the long term in order to clarify the impact of the
teleconnections. In figure 5, the evolution of the correlation skill between
the reference and the forecast of the variable $\psi_{a,1}$ with or without
errors in the initial conditions of the tropical model, is displayed. After a
quick decrease to the level of the PTP experiment (black continuous line), the
correlation skill persists around that level for a few days or weeks depending
on the amplitude of the error in the initial conditions of the tropical
forcing. It then drops to a lower level that also depends on the amplitude of
the initial error. Two important messages here: (i) The teleconnection skill
can never be reached on the long term basis when error in the initial
conditions are present in the tropical model forcing; (ii) the ability to keep
skill on long lead times depends on the amplitude of the errors in the initial
conditions of the tropical model.

Figure 5: Evolution of the correlation skill for $M=1$ for different
amplitudes of the errors introduced in the initial conditions of the tropical
model.

###  3.4 Impact of ensemble averaging

A first approach to improve the skill is to generate an ensemble of runs
starting from different initial conditions and to extract the ensemble mean.
It has been shown that the skill is improved even if a very small number of
ensemble members is used [38, 39].

Let us start by investigating the error evolution for different values of M
for the case for which the initial condition error amplitude of the tropical
model is fixed to $PT=0.05$. The results are displayed in figure (6). The key
information that can be extracted is that even with 5-10 members, the main
improvement that can be obtained through averaging is reached: little
improvements can be obtained going beyond that order. This was already pointed
out in [38] in an intermediate complexity model with much more degrees of
freedom.

Figure 6: MSE evolution of the forecast of $\psi_{a,1}$ for the control run
with $M=1$ (black continuous curve referred as control), and of the ensemble
mean for different values of the ensemble size from 5 to 50.

This impact of ensemble averaging on the correlation skill is then illustrated
in figure 7 for (a) $\psi_{a,1}$ with a perturbation amplitude of $PT=0.02$,
and (b) with a perturbation amplitude of $PT=0.05$. In these plots, the
correlation skill corresponding to the PTP is also displayed, as well as the
case with $M=1$.

A first remark is that the correlation skill between the reference and the
ensemble mean forecasts is still lower than PTP for long lead times when
initial condition perturbations are introduced in the forecasting model of the
forcing. Moreover when the number of members is increased, the quality of the
forecast is increased, but already 10 members are enough to reach the optimal
skill. The latter result is consistent with the fact that the ensemble mean is
rapidly converging to its asymptotic value as discussed in [38].

A second important result is the improvement of the correlation skill when the
ensemble mean is used as forecast between 50 and about 300 days. The
correlation skill is now reaching values close to the teleconnection
correlation between the forcing and the variable. In other words, the ensemble
averaging is properly filtering the fast weather variability around the mean.
The ensemble mean is here an appropriate approach to isolate the signal
present in the variable $\psi_{a,1}$.

Figure 7: Correlation skill of the forecasts for (a) $\psi_{a,1}$ with a
amplitude perturbation of $PT=0.02$, (b) $\psi_{a,1}$ with a amplitude
perturbation of $PT=0.05$. Within each panel, the curves are representing
different number of members M in the estimation of the ensemble mean. For
reference, the case without errors in the initial conditions of the tropical
model for $M=1$ is also displayed (black continuous curve).

For the variable $\theta_{a,1}$, the skill decreases quickly and is similar
whatever the initial condition error in the tropical model as illustrated in
figure 8, suggesting the weak sensitivity of this variable to the external
forcing. Note however that there is some improvements by taking ensemble
averages in the phase of quick decrease of the skill around 30-40 days, even
with a small number of ensemble members.

Figure 8: As in figure 7 but for the variable $\theta_{a,1}$ with an amplitude
perturbation of $PT=0.05$. The curves are representing different number of
members $M$ in the estimation of the ensemble mean. For reference the case
without errors in the initial conditions of the tropical model is also
displayed (black continuous curve).

In the last decades, many studies have tried to evaluate the potential
forecasting skill of climate forecasting models, and compared it to the actual
forecast skill [40, 41, 42]. The potential forecasting skill is the
correlation between the ensemble mean and any member of the ensemble. In this
context if the model is perfect and if the model is properly initialized, the
potential predictability should be equal to the actual predictability.

In fact even if the model is properly representing the reference evolution, a
discrepancy between the actual skill and the potential skill is found in the
current experiment. This is illustrated in figure 9. The potential skill is
larger than the actual skill, except when there is no initial condition errors
in the tropical model forecasts (black curves). This feature in the current
coupled tropical-extratropical model may only be understood by the fact that
even if perturbing randomly and symmetrically the initial conditions around
the observations (and hence the reference) in the tropical model, a
forecasting bias is emerging (not shown). This forecasting bias is most
probably due to the way the solutions of the forecasting model are converging
back to the model’s attractor. Such feature is very subtle and is worth
addressing further in a future studies by analyzing the temporal properties of
these transient trajectories, and to propose new approaches for perturbing the
model more consistent with the natural dynamics of the reference system as
discussed for instance in details in [43] or using analog initialization as in
[44].

Figure 9: Comparison of the actual forecasting skill of the ensemble mean and
the potential predictability – as defined by the skill of the ensemble mean in
forecasting any of the member of the ensemble. Three cases are displayed with
or without error in the initial conditions of the tropical model.

###  3.5 Impact of temporal averaging

Ensemble averages are often used to improve the estimate the the future state
of the system. As shown in the previous section, this is indeed interesting
and the correlation skill score can be improved to a certain extent. But once
the number of members is large enough, no improvement is obtained, this number
being related to the quick convergence of the ensemble mean to its asymptotic
value.

Another way to extract the signal from the forecasts is often to take temporal
averages (running means). This question has been addressed in several papers
in the past, in particular showing that the short term dynamics of the error
as described by the dominant Lyapunov exponent is the same as for
instantaneous fields [45, 46, 47]. The main modification that can then expect
is for long term forecasts as discussed in [48], but the gains are very low
when investigating atmosphere-only models [49]. We also explored this question
in the context of the tropical-extratropical coupled model to clarify the
impact on the time averaging in the extraction of the signal associated with
the tropical forcing.

Figure 10: MSE evolution for $\psi_{a,1}$ with $M=5$ and different temporal
averaging from $T=5.61$ d to $T=67.2$ d. Note that the lead time Time$=0$
corresponds to the end of the averaging period of the forecast.

Figure 10 is displaying the different MSE curves for the ensemble mean with
$M=5$ for the temporal averages of the variable $\psi_{a,1}$. The averages
have been performed from about 5 days to 60 days. Note that the lead time,
Time$=0$ days, is taken at the end of the averaging period $T$, so that
Time$=0$ on the figure is associated with the averaging over $T$ days of the
original forecast. Four different phases can be found in the evolution of the
error: first an overall increase of the MSE for short times associated with
the averaging of successive situations displaying large errors as compared
with the reference for longer lead times; a second phase during which the MSE
amplitude is decreasing when $T$ becomes large (from lead times of 25 to 200
days); a third phase of rapid growth (from 200 days to 400 days) for which
averaging does not change the MSE amplitude; and finally a phase (after 400
days) for which the MSE amplitude again decreases for large values of $T$.
This complicated dependence reflects that time averaging may indeed reduce the
variance, and hence the error, of the processes but not uniformly for all the
time scales. This aspect was already amply discussed in [45], and shows up
here in different responses of the MSE phases to averaging. This differential
behavior as a function of the phase does not however shows up for
$\theta_{a,1}$ as illustrated in figure 11 with a constant decrease of the MSE
saturation level as a function of time averaging. As $\theta_{a,1}$ is much
less affected by the forcing than $\psi_{a,1}$, this differential behavior for
$\psi_{a,1}$ could be associated with the action of the tropical forcing.

Figure 11: MSE evolution for $\theta_{a,1}$ with $M=5$ and different temporal
averaging from $T=5.61$ d to $T=67.2$ d. Note that the lead time $=0$
corresponds to the end of the averaging period of the forecast.

Another very useful view of what is happening while taking temporal averaging
is to evaluate the changes in the correlation skill. In figure 12, the
evolution of the correlation skill for the ensemble mean of 5 members ($M=5$)
of $\psi_{a,1}$ without time averaging, and with three different time averages
of 11.2, 33 and 67.2 days. A small gain of averaging is found between 50 and
300-day lead times, suggesting the usefulness of this averaging at such lead
times. If now one looks at the corresponding figure for $\theta_{a,1}$ (Fig.
13), the picture is very different. Although the instantaneous quantities,
with or without ensemble averaging, do not show skill beyond 40 days, the
temporal average allows for getting skill of this observable at long lead
times: For the temporal averaging of 67.2 days, a correlation skill of the
order of 0.2 is found up to almost 1000 days. This suggests that the temporal
averaging is appropriately removing the high frequency variability of
$\theta_{a,1}$ and isolate the signal associated with the influence of the
tropical forcing.

Figure 12: Evolution of the corrrelation skill for $\psi_{a,1}$ for various
combinations of averaging (ensemble and/or in time). Figure 13: Evolution of
the corrrelation skill for $\theta_{a,1}$ for various combinations of
averaging (ensemble and/or in time).

To sum up, temporal averages can be very helpful in extracting the skill of
observables influenced by external forcing signals, even when the
teleconnections are weak.

##  4 Conclusions

Low-frequency variability coming from the presence of an external slowly-
evolving forcing are often isolated through teleconnections. In the present
work, we have investigated in an idealized context the impact of the presence
of teleconnections between the tropics and the extratropics on the limits of
predictability of the extratropical fields. This has been done in a coupled
tropical-extratropical model that was recently built to evaluate the emergence
of LFV in the extratropics and the properties of the corresponding pullback
attractors [12].

A highly idealized twin experiment was performed by running the coupled
tropical-extratropical model with initial condition errors for both components
of the system, and comparing its forecasts to a long reference run. A few
generic properties were found:

  * •

If there is no error in the initial conditions of the tropical model, the
forced variable, here the first barotropic atmospheric mode, shows a skill on
very long time scales thanks to the slowly-varying external forcing (which is
in this case the same for the model forecast and the reference).

  * •

When errors in the initial conditions are present, the correlation skill is
considerably degraded at interannual lead times. The degradation of skill also
depends on the amplitude of the initial condition errors in the tropics,
suggesting that accurate initial conditions are necessary to get the full
potential of teleconnections.

  * •

Ensemble averaging may improve the quality of the forecasts, but the number of
members necessary to reach an appropriate skill does not need to be large:
5-10 members are enough in this perfect model framework.

  * •

Temporal averaging can also improve the quality of the forecasts (at the
expense of course of being able to forecast extreme events), in particular
when the forcing affects weakly the observable under interest. In the current
context, the first baroclinic mode shows a very weak teleconnection
(correlation) with the tropical forcing, and when averaging, a strong
correlation skill is emerging.

This work also showed that even if the model is perfect and the initial
observation distribution is known, the actual correlation skill of a forecast
is always less or equal than the potential correlation skill as defined using
the forecasting model itself only. This feature is very subtle and is most
probably associated with the way the model forecasts are converging back to
the system’s attractor. This attractor is the same as the reference, but the
path to converge to it is generating a forecasting bias which reduces the
actual correlation skill. This suggests that appropriate initial conditions of
the tropical model generating the forcing should be looked for, allowing for
reducing this bias. This will be the subject of a future work along the lines
of [43].

In the current work, the low-frequency variability is introduced through a one
way external forcing from the tropics to the midlatitudes. Things become more
complicated when part of the low-frequency variability is generated
intrinsically in the forced system as for instance in [50] or when the forcing
is interacting in a complicated way with the intrinsic low-frequency
variability like in [12]. In both cases, the mean may loose significance, and
a full probability distribution is then necessary to characterize the dynamics
of the reference system. Further investigations of such complex interactions
are worth performing in the future.

Finally this problem should also be tackled in the context of a more realistic
model for which the initial conditions in the tropical domain should be
modified in order to clarify how well we can exploit the teleconnections for
prediction purposes. This will also be the subject of a future work.

This work is supported by the Belgian Policy Office under Grant
B2/20E/P1/ROADMAP, a contribution to the ROADMAP project of the JPI-Climate-
JPI-Ocean research activities.

## References

  * [1] S.G.H. Philander.  El Niño and the Southern Oscillation.  Academic Press, 1990.
  * [2] M. A. Alexander, I. Bladé, M. Newman, J. R. Lanzante, N-Ch. Lau, and J. D. Scott.  The atmospheric bridge: The influence of ENSO teleconnections on air–sea interaction over the global oceans.  J. Climate, 15:2205–2231, 2002.
  * [3] M. P. Hoerling and A. Kumar.  Atmospheric response patterns associated with tropical forcing.  J. Climate, 15:2184–2203, 2002.
  * [4] C. Stan, D. M. Straus, J. S. Frederiksen, H. Lin, E. D. Maloney, and C. Schumacher.  Review of tropical-extratropical teleconnections on intraseasonal time scales.  Reviews of Geophysics, 55(4):902–937, 2017.
  * [5] S-W. Yeh, W. Cai, S-K. Min, M. J. McPhaden, D. Dommenget, B. Dewitte, M. Collins, K. Ashok, S-I. An, B-Y. Yim, and J-S. Kug.  Enso atmospheric teleconnections and their response to greenhouse gas forcing.  Reviews of Geophysics, 56(1):185–206, 2018.
  * [6] T. J. Mosedale, D. B. Stephenson, M. Collins, and T. C. Mills.  Granger causality of coupled climate processes: Ocean feedback on the north atlantic oscillation.  Journal of Climate, 19(7):1182 – 1194, 2006.
  * [7] M. Paluš, A. Krakovská, J. Jakubík, and M. Chvosteková.  Causality, dynamical systems and the arrow of time.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 28(7):075307, 2018.
  * [8] J. Runge, S. Bathiany, E. Bollt, G. Camps-Valls, D. Coumou, E. Deyle, C. Glymour, M. Kretschmer, M. D. Mahecha, J. Muñoz-Marí, E. H. van Nes, J. Peters, R. Quax, M. Reichstein, M. Scheffer, B. Schölkopf, P. Spirtes, G. Sugihara, J. Sun, K. Zhang, and J. Zscheischler.  Inferring causation from time series in earth system sciences.  Nature Communications, 10(1):2553, 2019.
  * [9] S. Vannitsem and X. S. Liang.  Dynamical dependencies at monthly and interannual time scales in the climate system: Study of the north pacific and atlantic regions.  Tellus A: Dynamic Meteorology and Oceanography, Mar 2022.
  * [10] A. Kumar and M. P. Hoerling.  Prospects and limitations of seasonal atmospheric gcm predictions.  Bull. Amer. Meteor. Soc., 76:335–345, 1995.
  * [11] A. G. Nidheesh, M. Lengaigne, J. Vialard, T. Izumo, A. S. Unnikrishnan, and Ch. Cassou.  Influence of ENSO on the Pacific decadal oscillation in CMIP models.  Climate Dynamics, 49:3309–3326, 2017.
  * [12] S. Vannitsem, J. Demaeyer, and M. Ghil.  Extratropical low-frequency variability with enso forcing: A reduced-order coupled model study.  Journal of Advances in Modeling Earth Systems, 13(6):e2021MS002530, 2021.  e2021MS002530 2021MS002530.
  * [13] A. Timmermann, F-F. Jin, and J. Abshagen.  A nonlinear theory for El Niño bursting.  J. Atmos. Sci., 60:152–165, 2003.
  * [14] S. Vannitsem.  The role of the ocean mixed layer on the development of the north atlantic oscillation: A dynamical system’s perspective.  Geophysical Research Letters, 42(20):8615–8623, 2015.
  * [15] F-F. Jin.  Tropical ocean-atmosphere interaction, the Pacific cold tongue, and the El-Niño-Southern Oscillation.  Science, 274:76–78, 1996.
  * [16] F-F. Jin.  An equatorial ocean recharge paradigm for ENSO. Part I: Conceptual model.  J. Atmos. Sci., 54:811–829, 1997.
  * [17] S-I. An and F-F. Jin.  Nonlinearity and asymmetry of ENSO.  J. Climate, 17:2399–2412, 2004.
  * [18] A. Roberts, J. Guckenheimer, E. Widiasih, A. Timmermann, and Ch. K. R. T. Jones.  Mixed-mode oscillations of El Niño-Southern Oscillation.  J Atmos Sci, 73:1755–1766, 2016.
  * [19] S. Vannitsem, J. Demaeyer, L. De Cruz, and M. Ghil.  Low-frequency variability and heat transport in a low-order nonlinear coupled ocean-atmosphere model.  Physica D: Nonlinear Phenomena, 309:71–85, 2015.
  * [20] L. De Cruz, J. Demaeyer, and S. Vannitsem.  The Modular Arbitrary-Order Ocean-Atmosphere Model: MAOOAM v1.0.  Geosci. Mod. Dev., 9(8):2793–2808, 2016.
  * [21] S. Vannitsem.  Predictability of large-scale atmospheric motions: Lyapunov exponents and error dynamics.  Chaos, 27:32101, 2017.
  * [22] S. Vannitsem and V. Lucarini.  Statistical and dynamical properties of covariant Lyapunov vectors in a coupled atmosphere-ocean model – multiscale effects, geometric degeneracy, and error dynamics.  J. Phys. A, 49(22):224001, 2016.
  * [23] L. De Cruz, S. Schubert, J. Demaeyer, V. Lucarini, and S. Vannitsem.  Exploring the Lyapunov instability properties of high-dimensional atmospheric and climate models.  Nonlin. Proc. Geophys., 25(2):387–412, 2018.
  * [24] J. Demaeyer and S. Vannitsem.  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(703):881–896, 2017.
  * [25] J. Demaeyer and S. Vannitsem.  Comparison of stochastic parameterizations in the framework of a coupled ocean-atmosphere model.  Nonl. Proc. Geophys., 25(3):605–631, 2018.
  * [26] S.G. Penny, E. Bach, K. Bhargava, C.-C. Chang, C. Da, L. Sun, and T. Yoshida.  Strongly coupled data assimilation in multiscale media: Experiments using a quasi-geostrophic coupled model.  J. Adv. Model. Earth Syst, 6:1803–1829, 2019.
  * [27] M. Tondeur, A. Carrassi, S. Vannitsem, and M. Bocquet.  On temporal scale separation in coupled data assimilation with the ensemble Kalman filter.  J Stat Phys, 179:1161–1185, 2020.
  * [28] A. Carrassi, C. Grudzien, M. Bocquet, J. Demaeyer, P. Raanes, and S. Vannitsem.  Data assimilation for chaotic systems.  In Seon-Ki Park and Xu Liang, editors, Data Assimilation for Atmospheric, Oceanic and Hydrological Applications. Springer Science & Business Media, 2022.
  * [29] A. E. Gill.  Atmosphere-Ocean Dynamics.  Academic Press, New York, U.S.A., 1982.
  * [30] J. Pedlosky.  Geophysical Fluid Dynamics.  Springer-Verlag, New York, 2nd edition, 1987.
  * [31] S. Schemm, G. Rivière, L. M. Ciasto, and C. Li.  Extratropical cyclogenesis changes in connection with tropospheric ENSO teleconnections to the North Atlantic: Role of stationary and transient waves.  J. Atmos. Sci., 75:3943–3964, 2018.
  * [32] Andrew T. Wittenberg.  Are historical records sufficient to constrain enso simulations?  Geophysical Research Letters, 36(12), 2009.
  * [33] John Guckenheimer, Axel Timmermann, Henk Dijkstra, and Andrew Roberts.  (Un)predictability of strong El Niño events.  Dynamics and Statistics of the Climate System, 2(1):dzx004, 12 2017\.
  * [34] M. Ghil, M. R. Allen, M. D. Dettinger, K. Ide, D. Kondrashov, M. E. Mann, A. W. Robertson, A. Saunders, Y. Tian, F. Varadi, and P. Yiou.  Advanced spectral methods for climatic time series.  Reviews of Geophysics, 40(1):41 pages, 2002.
  * [35] S. Vannitsem and Z. Toth.  Short-term dynamics of model errors.  Journal of the Atmospheric Sciences, 59(17):2594 – 2604, 2002.
  * [36] C. Nicolis.  Dynamics of model error: Some generic features.  Journal of the Atmospheric Sciences, 60(17):2208 – 2218, 2003.
  * [37] C. Nicolis.  Dynamics of model error: The role of the boundary conditions.  Journal of the Atmospheric Sciences, 64(1):204 – 215, 2007.
  * [38] P. L. Houtekamer and J. Derome.  Methods for ensemble prediction.  Monthly Weather Review, 123(7):2181 – 2196, 1995.
  * [39] M. Ehrendorfer.  Predicting the uncertainty of numerical weather forecasts: A review.  METEOROLOGISCHE ZEITSCHRIFT-BERLIN-, 6:147–183, 1997.
  * [40] A. Kumar, P. Peng, and M. Chen.  Is there a relationship between potential and actual skill?  Monthly Weather Review, 142(6):2220 – 2227, 2014.
  * [41] Y. Imada and H. Kawase.  Potential seasonal predictability of the risk of local rainfall extremes estimated using high-resolution large ensemble simulations.  Geophysical Research Letters, 48(24):e2021GL096236, 2021.  e2021GL096236 2021GL096236.
  * [42] D. Yang, Y. Tang, XQ. Yang, X. Song, X. Tan, Y. Wu, X. Yan, T. Liu, and X. Sun.  Probabilistic versus deterministic potential seasonal climate predictability under the perfect-model framework.  Climate Dynamics, 2023.
  * [43] J. Demaeyer, S. G. Penny, and S. Vannitsem.  Identifying efficient ensemble perturbations for initializing subseasonal-to-seasonal prediction.  Journal of Advances in Modeling Earth Systems, 14(5):e2021MS002828, 2022.  e2021MS002828 2021MS002828.
  * [44] S. Li, X. Rong, Y. Liu, Z. Liu, and K. Fraedrich.  Dynamic analogue initialization for ensemble forecasting.  Advances in Atmospheric Sciences, 30(5):1406–1420, 2013.
  * [45] C. Nicolis and G. Nicolis.  From short-scale atmospheric variability to global climate dynamics: Toward a systematic theory of averaging.  Journal of Atmospheric Sciences, 52(11):1903 – 1913, 1995.
  * [46] S. Vannitsem and C. Nicolis.  Lyapunov vectors and error growth patterns in a t21l3 quasigeostrophic model.  Journal of the Atmospheric Sciences, 54(2):347 – 361, 1997.
  * [47] S. Vannitsem and C. Nicolis.  Dynamics of fine-scale variables versus averaged observables in a t21l3 quasi-geostrophic model.  Quarterly Journal of the Royal Meteorological Society, 124(551):2201–2226, 1998.
  * [48] J.O. Roads.  Predictability in the extended range.  Journal of Atmospheric Sciences, 44(23):3495 – 3527, 1987.
  * [49] J. J. Tribbia and D. P. Baumhefner.  Estimates of the predictability of low-frequency variability with a spectral general circulation model.  Journal of Atmospheric Sciences, 45(16):2306 – 2318, 1988.
  * [50] S. Pierini.  Statistical significance of small ensembles of simulations and detection of the internal climate variability: An excitable ocean system case study.  J. Stat. Phys., 179:1475–1495, 2020.