# The predictability problem in systems with an uncertainty in the evolution law

## I Introduction

The ability to predict the future state of a system, given its present state, stands at the foundations of scientific knowledge with relevant implications from an applicative point of view in geophysical and astronomical sciences. In the prediction of the evolution of a system, e.g. the atmosphere, we are severely limited by the fact that we do not know with arbitrary accuracy the evolution equations and the initial conditions of the system. Indeed, one integrates a mathematical model given by a finite number of equations. The initial condition, a point in the phase space of the model, is determined only with a finite resolution (i.e. by a finite number of observations) (Monin 1973).

Using concepts of dynamical systems theory, there have been some progresses in understanding the growth of an uncertainty during the time evolution. An infinitesimal initial uncertainty ($`\delta _00`$) in the limit of long times ($`t\mathrm{}`$) grows exponentially in time with a typical rate given by the leading Lyapunov exponent $`\lambda `$, $`|\delta x(t)|\delta _0\mathrm{exp}(\lambda t)`$. Therefore if our purpose is to forecast the system within a tolerance $`\mathrm{\Delta }`$, the future state of the system can be predicted only up to the predictability time, given by:

$$T_p\frac{1}{\lambda }\mathrm{ln}\left(\frac{\mathrm{\Delta }}{\delta _0}\right).$$
(1)   
In literature, the problem of predictability with respect to uncertainty on the initial conditions is referred to as predictability of the first kind.

In addition, in real systems we must also cope with the lack of knowledge of the evolution equations. Let us consider a system described by a differential equation:

$$\frac{d}{dt}𝐱(t)=𝐟(𝐱,t),𝐱,𝐟^n.$$
(2)   
As a matter of fact we do not know exactly the equations, and we have to devise a model which is different from the true dynamics:

$$\frac{d}{dt}𝐱(t)=𝐟_ϵ(𝐱,t)\text{where}𝐟_ϵ(𝐱,t)=𝐟(𝐱,t)+ϵ\delta 𝐟(𝐱,t).$$
(3)   
Therefore, it is natural to wonder about the relation between the true evolution (reference or true trajectory $`𝐱_T(t)`$) given by (2) and that one effectively computed (perturbed or model trajectory $`𝐱_M(t)`$) given by (3). This problem is referred to as predictability of the second kind.

Let us make some general remarks. At the foundation of the second kind predictability problem there is the issue of structural stability (Guckenheimer et al. 1983): since the evolution laws are known only with finite precision it is highly desirable that at least certain properties are not too sensitive to the details of the equations of motion. For example, in a system with a strange attractor, small generic changes in the evolution laws should not change drastically the dynamics (see Appendix A for a simple example with non generic perturbation).

In chaotic systems the effects of a small generic uncertainty on the evolution law are similar to those due to the finite precision on the initial condition (Crisanti et al. 1989). The model trajectory of the perturbed dynamics diverges exponentially from the reference one with a mean rate given by the Lyapunov exponent of the original system. The statistical properties (such as correlation functions and temporal averages) are not strongly modified. This last feature has been frequently related to the shadowing lemma (Guckenheimer et al. 1983; Ott 1993): almost all trajectories of the true system can be approximated by a trajectory of the perturbed system starting from a slightly different initial condition. However, as far as we know, the shadowing lemma can be proven only in special cases and therefore it cannot be straightforwardly invoked to explain the statistical reproducibility in a generic case. In addition, in real systems the size of an uncertainty on the evolution equations is determinable only a posteriori, based on the ability of the model equations to reproduce some of the features of the phenomenon.

In dynamical systems theory, the problems of first and second kind predictability is essentially understood in the limit of infinitesimal perturbations. However even in this limit we must also consider the fluctuations of the rate of expansion which can lead to relevant modifications of the predictability time (1), in particular for strongly intermittent systems (Benzi et al. 1985; Paladin et al. 1987; Crisanti et al. 1993).

As far as finite perturbations are considered, the leading Lyapunov exponent is not relevant for the predictability issue. In presence of many characteristic times and spatial scales the Lyapunov exponent is related to the growth of small scale perturbations which saturates on short times and has very little relevance for the growth of large scale perturbations (Leith and Kraichnan 1972; Monin 1973; Lorenz 1996). To overcome this shortcoming, a suitable characterization of the growth of non infinitesimal perturbations, in terms of the Finite Size Lyapunov Exponent (FSLE), has been recently introduced (Aurell et al. 1996 and 1997).

Also in the case of second kind predictability one has often to deal with errors which are far from being infinitesimal. Typical examples are systems described by partial differential equations (e.g. turbulence, atmospheric flows). The study of these systems is performed by using a numerical model with unavoidable severe approximations, the most relevant of which is the necessity to cut some degrees of freedom off; basically, the small scale variables.

The aim of this paper is to analyze the effects of limited resolution on the large scale features. This raises two problems: in first place one has to deal with perturbations of the evolution equations which in general cannot be considered small; second, the parameterization of the unresolved modes can be a subtle point. We shall show that the Finite Size Lyapunov Exponent is able to characterize the effects of uncertainty on the evolution laws. Moreover we shall discuss the typical difficulties arising in the parameterization of the unresolved scales.

This paper is organized as follows. In section II we report some known results about the predictability problem of the second kind and recall the definition of the FSLE. In section III we present numerical results on simple models. In section IV we consider more complex systems with many characteristic times. Section V is devoted to summarize the results. In Appendix A we illustrate a simple example of structural unstable system. In Appendix B we describe the method for the computation of the FSLE and in Appendix C we discuss the problem of the parameterization of the unresolved variables.

## II EFFECTS OF A SMALL UNCERTAINTY ON THE EVOLUTION LAW

In the second kind predictability problem, we can distinguish three general cases depending on the original dynamics. In particular, equation (2) may display:

(i) trivial attractors: asymptotically stable fixed points or attracting periodic orbits;

(ii) marginally stable fixed points or periodic/quasi-periodic orbits as in integrable Hamiltonian systems;

(iii) chaotic behavior.

In case (i) small changes in the equations of motion do not modify the qualitative features of the dynamics. Case (ii) is not generic and the outcome strongly depends on the specific perturbation $`\delta 𝐟`$, i.e. it is not structurally stable. In the chaotic case (iii) one expects that the perturbed dynamics is still chaotic. In this paper we will consider only this latter case.

Let us also mention that, in numerical computations of evolution equations (e.g. differential equations), there are two unavoidable sources of errors: the finite precision representation of the numbers which causes the computer phase space to be necessarily discrete and the round-off which introduces a sort of noise. Because of the discrete nature of the phase space of the system studied on computer, orbits numerically computed have to be periodic. Nevertheless the period is usually very large, apart for very low computer precision (Crisanti et. al 1989). We do not consider here this source of difficulties. The round-off produces on eq. (2) a perturbation which can be written as $`\delta 𝐟(𝐱,t)=𝐰(𝐱)𝐟(𝐱,t)`$ and $`ϵ10^\alpha `$ ($`\alpha =`$number of digits in floating point representation) where $`𝐰=O(1)`$ is an unknown function which may depend on $`𝐟`$ and on the software of the computer (Knut 1969). In general, the round-off error is very small and may have, as much as the noise, a positive role, as underlined by Ruelle (1979), in selecting the physical probability measure, the so-called natural measure, from the set of ergodic invariant measures.

In chaotic systems the effects of a small uncertainty on the evolution law is, for many aspects, similar to those due to imperfect knowledge of initial conditions. This can be understood by the following example. Consider the Lorenz equations (Lorenz 1963)

$$\begin{array}{ccc}\frac{dx}{dt}\hfill & =\hfill & \sigma (yx)\hfill \\ \frac{dy}{dt}\hfill & =\hfill & Rxyxz\hfill \\ \frac{dz}{dt}\hfill & =\hfill & xybz.\hfill \end{array}$$
(4)   
In order to mimic an experimental error in the determination of the evolution law we consider a small error $`ϵ`$ on the parameter $`R`$: $`RR+ϵ`$. Let us consider the difference $`\mathrm{\Delta }𝐱(t)=𝐱_M(t)𝐱_T(t)`$ with, for simplicity, $`\mathrm{\Delta }𝐱(0)=0`$, i.e. we assume a perfect knowledge of the initial conditions. One has, with obvious notation:

$$\frac{d\mathrm{\Delta }𝐱}{dt}=𝐟_ϵ(𝐱_M)𝐟(𝐱_T)\frac{𝐟}{𝐱}\mathrm{\Delta }𝐱+\frac{𝐟_ϵ}{R}ϵ.$$
(5)   
At time $`t=0`$ one has $`|\mathrm{\Delta }𝐱(0)|=0`$, therefore $`|\mathrm{\Delta }𝐱(t)|`$ grows initially only by the effect of the second term in (5). At later times, when $`|\mathrm{\Delta }𝐱(t)|O(ϵ)`$ the leading term of (5) becomes the first one, and we recover the first kind predictability problem for an initial uncertainty $`\delta _0ϵ`$. Therefore, apart from an initial (not particularly interesting) growth, which depends strongly on the specific perturbation, the evolution of $`<\mathrm{log}(|\mathrm{\Delta }𝐱(t)|)>`$ follows the usual linear growth with the slope given by the leading Lyapunov exponent. Typically the value of the Lyapunov exponent computed by using the model dynamics differs from the true one by a small amount of order $`ϵ`$, i.e. $`\lambda _M=\lambda _T+O(ϵ)`$ (Crisanti et al. 1989).

This consideration applies only to infinitesimal perturbations. The generalization to finite perturbations requires the extension of the Lyapunov exponent to finite errors. Let us now introduce the Finite Size Lyapunov Exponent for the predictability of finite perturbations. The definition of FSLE $`\lambda (\delta )`$ is given in terms of the “doubling time” $`T_r(\delta )`$, that is the time a perturbation of initial size $`\delta `$ takes to grow by a factor $`r`$ ($`>1`$):

$$\lambda (\delta )=\frac{1}{T_r(\delta )}_t\mathrm{ln}r$$
(6)   
where $`\mathrm{}_t`$ denotes average with respect to the natural measure, i.e. along the trajectory (see Appendix B). For chaotic systems, in the limit of infinitesimal perturbations ($`\delta 0`$) $`\lambda (\delta )`$ is nothing but the leading Lyapunov exponent $`\lambda `$ (Benettin et al. 1980). Let us note that the above definition of $`\lambda (\delta )`$ is not appropriate to discriminate cases with $`\lambda =0`$ and $`\lambda <0`$, since the predictability time is positive by definition. Nevertheless this is not a limitation as long as we deal with chaotic systems.

In many realistic situations the error growth for infinitesimal perturbations is dominated by the fastest scales, which are typically the smallest ones (e.g. small scale turbulence). When $`\delta `$ is no longer infinitesimal, $`\lambda (\delta )`$ is given by the fully nonlinear evolution of the perturbation. In general $`\lambda (\delta )\lambda `$, according to the intuitive picture that large scales are more predictable. Outside the range of scales in which the error $`\delta `$ can be considered infinitesimal, the function $`\lambda (\delta )`$ depends on the details of the dynamics and in principle on the norm used. In fully developed turbulence one has the universal law $`\lambda (\delta )\delta ^2`$ in the inertial range (Aurell et al. 1996 and 1997). It is remarkable that this prediction, which can be obtained within the multifractal model for turbulence, is not affected by intermittency and it gives the law originally proposed by Lorenz (1969). The behavior of $`\lambda (\delta )`$ as a function of $`\delta `$ gives important information on the characteristic times and scales of the system and it has been also applied to passive transport in closed basins (Artale et al. 1997).

Let us now return to the example (4). We compute $`\lambda _{TT}(\delta )`$, the FSLE for the true equations, and $`\lambda _{TM}(\delta )`$, the FSLE computed following the distance between one true trajectory and one model trajectory starting at the same point. These are shown in Figure 2. The true FSLE $`\lambda _{TT}(\delta )`$ displays a plateau indicating a chaotic dynamics with leading Lyapunov exponent $`\lambda 1`$. Concerning the second kind predictability, for $`\delta >ϵ`$ the second term in (5) becomes negligible and we observe the transition to the Lyapunov exponent $`\lambda _{TM}(\delta )\lambda _{TT}(\delta )\lambda `$. In this range of errors the model system recovers the intrinsic predictability of the true system. For very small errors, $`\delta <ϵ`$, where the growth of the error is dominated by the second term in (5), we have $`\lambda _{TM}(\delta )>\lambda _{TT}(\delta )`$.

This example shows that it is possible to recover the intrinsic predictability of a chaotic system even in presence of some uncertainty in the model equations.

The relevance of the above example is however limited by the fact that (4) does not involve different scales. To investigate the effect of spatial resolution on predictability let us consider the advection of Lagrangian tracers in a given Eulerian field. We study a time-dependent, two dimensional, velocity field given by the superposition of large scale (resolved) eddies and small scale (possibly unresolved) eddies.

The streamfunction we consider is a slight modification of a model originally proposed for chaotic advection in Rayleigh-Bénard convection (Solomon et al. 1988):

$$\mathrm{\Psi }(x,y,t)=\psi (x,y,t;k_L,\omega _L,B_L)+ϵ\psi (x,y,t;k_S,\omega _S,B_S)$$
(7)   
with

$$\psi (x,y,t;k,\omega ,B)=\frac{1}{k}\mathrm{sin}\left\{k\left[x+B\mathrm{sin}(\omega t)\right]\right\}\mathrm{sin}\left(ky\right)$$
(8)   
The first term represents the large-scale flow, i.e. the resolved part of the flow, the second one mimics the unresolved small scale term and $`ϵ`$ measures the relative amplitude. We choose $`k_Sk_L`$ and $`\omega _S\omega _L`$ in order to have a sharp separation of space and time scales.

The Lagrangian tracers evolve according to the equations:

$$\frac{dx}{dt}=\frac{\psi }{y},\frac{dy}{dt}=\frac{\psi }{x}.$$
(9)   
We use the complete stream function (7) for the true dynamics and only the large-scale term for the model dynamics.

The time dependence induces chaotic motion and diffusion in the $`x`$ direction, without inserting any noise term (Solomon et al. 1988). For what concerns $`\lambda _{TT}`$ one observes three regimes (Figure 3). For very small errors $`\delta <2\pi /k_S`$ the exponential separation is ruled by the fastest scale and $`\lambda _{TT}\lambda `$, i.e. we recover the Lyapunov exponent of the system. At intermediate errors we observe a second small plateau corresponding to the large-scale term for the model dynamics. For larger errors, $`\delta >2\pi /k_L`$ one has $`\lambda _{TT}\delta ^2`$, i.e. diffusive behavior (see Artale et. al 1997).

The model FSLE $`\lambda _{TM}(\delta )`$ cannot recover the small scale features: for $`\delta ϵ`$ we observe the scaling $`\lambda _{TM}(\delta )\delta ^1`$ which can be understood by the following argument. In this region the distance between the reference and the true trajectories grows as $`d\delta /dtϵ`$ and thus, by a dimensional estimate, one has:

$$\lambda _{TM}(\delta )\frac{1}{T_r(\delta )}\frac{ϵ}{\delta }.$$
(10)   
Nevertheless, for larger $`\delta `$ the model fairly captures the small plateau displayed by $`\lambda _{TT}(\delta )`$ which corresponds to the slow time scale, then at large $`\delta `$ (i.e. for $`\delta `$ greater than the $`\delta >2\pi /k_L`$) we recover the diffusive behavior with the correct diffusion coefficient. This last feature can be understood by the fact that the diffusion coefficient, being an asymptotic quantity of the flow, is not influenced by the details of small scale structures.

This example is rather simple: large scales do not interact with the small ones and the number of degrees of freedom is very small. Therefore in this case the crude elimination of the small scale component does not prevent the possibility of a fair description of large scale features.

In the following we will consider more complex situations, in which strongly interacting degrees of freedom with different characteristic times are involved. In these cases the correct parameterization of the unresolved modes is crucial for the prediction of large scale behavior.

## III Systems with two time scales

Before analyzing in detail the effects of non infinitesimal perturbations of the evolution laws in some specific models let us clarify our aims. We consider a dynamical system written in the following form:

$$\begin{array}{ccc}\frac{d𝐱}{dt}\hfill & =\hfill & 𝐟(𝐱,𝐲)\hfill \\ \frac{d𝐲}{dt}\hfill & =\hfill & 𝐠(𝐱,𝐲),\hfill \end{array}$$
(11)   
where $`𝐟,𝐱^n`$ and $`𝐠,𝐲^m`$, in general $`nm`$. Now, let us suppose that the fast variables $`𝐲`$ cannot be resolved: a typical example are the subgrid modes in PDE discretizations. In this framework, a natural question is: how must we parameterize the unresolved modes ($`𝐲`$) in order to predict the resolved modes ($`𝐱`$) ?

As discussed by Lorenz (1996), to reproduce – at a qualitative level – a given phenomenology, e.g. the ENSO phenomenon, one can drop out the small scale features without negative consequences. But one unavoidably fails in forecasting the ENSO (i.e. the actual trajectory) without taking into account in a suitable way the small scale contributions.

An example in which it is relatively simple to develop a model for the fast modes is represented by skew systems:

$$\begin{array}{ccc}\frac{d𝐱}{dt}\hfill & =\hfill & 𝐟(𝐱,𝐲)\hfill \\ \frac{d𝐲}{dt}\hfill & =\hfill & 𝐠(𝐲)\hfill \end{array}$$
(12)   
In this case, the fast modes $`(𝐲)`$ do not depend on the slow ones $`(𝐱)`$. One can expect that in this case, neglecting the fast variables or parameterizing them with a suitable stochastic process, should not drastically affect the prediction of the slow variables (Boffetta et al. 1996).

On the other hand, if $`𝐲`$ feels some feedback from $`𝐱`$, we cannot simply neglect the unresolved modes. In Appendix C we discuss this point in detail. In practice one has to construct an effective equation for the resolved variables:

$$\frac{d𝐱}{dt}=𝐟_M(𝐱,𝐲(𝐱)),$$
(13)   
where the functional form of $`𝐲(𝐱)`$ and $`𝐟_M`$ are found by phenomenological arguments and/or by numerical studies of the full dynamics.

Let us now investigate an example with a recently introduced toy model of the atmosphere circulation (Lorenz 1996; Lorenz et al. 1998) including large scales $`x_k`$ (synoptic scales) and small scales $`y_{j,k}`$ (convective scales):

$$\begin{array}{ccc}\frac{dx_k}{dt}\hfill & =\hfill & x_{k1}\left(x_{k2}x_{k+1}\right)\nu x_k+F_{j=1}^Jy_{j,k}\hfill \\ \frac{dy_{j,k}}{dt}\hfill & =\hfill & cby_{j+1,k}\left(y_{j+2,k}y_{j1,k}\right)c\nu y_{j,k}+x_k\hfill \end{array}$$
(14)   
where $`k=1,\mathrm{},K`$ and $`j=1,\mathrm{},J`$. As in (Lorenz 1996) we assume periodic boundary conditions on $`k`$ ($`x_{K+k}=x_k`$, $`y_{j,K+k}=y_{j,k}`$) while for $`j`$ we impose $`y_{J+j,k}=y_{j,k+1}`$. The variables $`x_k`$ represent some large scale atmospheric quantities in $`K`$ sectors extending on a latitude circle, while the $`y_{j,k}`$ represent quantities on smaller scales in $`JK`$ sectors. The parameter $`c`$ is the ratio between fast and slow characteristic times and $`b`$ measures the relative amplitude.

As pointed out by Lorenz, this model shares some basic properties with more realistic models of the atmosphere. In particular, the non-linear terms, which model the advection, are quadratic and conserve the total kinetic energy $`_k(x_k^2+_jy_{j,k}^2)`$ in the unforced ($`F=0`$), inviscid ($`\nu =0`$) limit; the linear terms containing $`\nu `$ mimic dissipation and the constant term $`F`$ acts as an external forcing preventing the total energy from decaying.

If one is interested in forecasting the large scale behavior of the atmosphere by using only the slow variables, a natural choice for the model equations is:

$$\frac{dx_k}{dt}=x_{k1}\left(x_{k2}x_{k+1}\right)\nu x_k+FG_k(𝐱),$$
(15)   
where $`G_k(𝐱)`$ represents the parameterization of the fast components in (14) (see Appendix C).

The FSLE for the true system (Boffetta et al. 1998) is shown in Figure 4 and displays the two characteristic plateau corresponding to fast component (for $`\delta 0.1`$) and slow component for large $`\delta `$ . Figure 4 also shows what happens when one simply neglects the fast components $`y_{j,k}`$ (i.e. $`𝐆(𝐱)=0`$). At very small $`\delta `$ one has $`\lambda _{TM}(\delta )\delta ^1`$ as previously discussed. For large errors we observe that, with this rough approximation, we are not able to capture the characteristic predictability of the original system. More refined parameterizations in terms of stochastic processes with the correct probability distribution function and correlation times do not improve the forecasting ability.

The reason for this failure is due to the presence of a feedback term in the equations (14) which induces strong correlations between the variable $`x_k`$ and the unresolved coupling $`_{j=1}^Jy_{j,k}`$. For a proper parameterization of the unresolved variables we follow the strategy discussed in Appendix C. Basically we adopt

$$G(𝐱)=\nu _ex_k,$$
(16)   
in which $`\nu _e`$ is a numerically determined parameter. Figure 4 shows that, although small scale are not resolved, the large scale predictability is well reproduced and one has $`\lambda _{TM}(\delta )\lambda _{TT}(\delta )`$ for large $`\delta `$. We conclude this section by observing that the proposed parameterization (16) is a sort of eddy viscosity parameterization.

## IV Large scale predictability in a turbulence model

We now consider a more complex system which mimics the energy cascade in fully developed turbulence. The model is in the class of the so called shell models introduced some years ago for a dynamical description of small-scale turbulence. For a recent review on shell models see Bohr et al. 1998. This model has relatively few degrees of freedom but involves many characteristic scales and times. The velocity field is assumed isotropic and it is decomposed on a finite set of complex velocity components $`u_n`$ representing the typical turbulent velocity fluctuation on a “shell” of scales $`\mathrm{}_n=1/k_n`$. In order to reach very high Reynolds number with a moderate number of degrees of freedom, the scales are geometrically spaced as $`k_n=k_02^n`$ ($`n=1,\mathrm{}N`$).

The specific model here considered has the form (L’vov et al. 1998)

$$\frac{du_n}{dt}=i\left(k_{n+1}u_{n+1}^{}u_{n+2}\frac{1}{2}k_nu_{n1}^{}u_{n+1}+\frac{1}{2}k_{n1}u_{n2}u_{n1}\right)\nu k_n^2u_n+f_n$$
(17)   
where $`\nu `$ represent the kinematic viscosity and $`f_n`$ is a forcing term which is restricted only to the first two shells (in order to mimic large scale energy injection).

Without entering in the details, we recall that the Shell Model (17) displays an energy cascade á la Kolmogorov from large scales (small $`n`$) to dissipative scales ($`nN`$) with a statistical stationary energy flux. Scaling laws for the average velocity components are observed:

$$|u_n^p|k_n^{\zeta _p}$$
(18)   
with exponents close to the Kolmogorov 1941 values $`\zeta _p=p/3`$.

From a dynamical point of view, model (17) displays complex chaotic behavior which is responsible of the small deviation of the scaling exponents (intermittency) with respect to the Kolmogorov values. Neglecting this (small) intermittency effects, a dimensional estimate of the characteristic time (eddy turnover time) for scale $`n`$ gives

$$\tau _n\frac{\mathrm{}_n}{|u_n|}k_n^{2/3}.$$
(19)   
The scaling behavior holds up to the Kolmogorov scale $`\eta =1/k_d`$ defined as the scale at which the dissipative term in (17) becomes relevant. The Lyapunov exponent of the turbulence model can be estimated as the fastest characteristic time $`\tau _d`$ and one has the prediction (Ruelle 1979)

$$\lambda \frac{1}{\tau _d}Re^{1/2}$$
(20)   
where we have introduced the Reynolds number $`Re1/\nu `$. It is possible to predict the behavior of the FSLE by observing that the faster scale $`k_n`$ at which an error of size $`\delta `$ is still active (i.e. below the saturation) is such that $`u_n\delta `$. Thus $`\lambda (\delta )1/\tau _n`$ and, using Kolmogorov scaling, one obtains (Aurell et al. 1996 and 1997)

$$\lambda _{TT}(\delta )\{\begin{array}{ccc}\lambda \hfill & \text{for}& \delta u_d\hfill \\ \delta ^2\hfill & \text{for}& u_d\delta u_0\hfill \end{array}$$
(21)   
To be more precise there is an intermediate range between the two showed in (21). For a discussion on this point see (Aurell et al. 1996 and 1997).

In order to simulate a finite resolution in the model, we consider a modelization of (17) in terms of an eddy viscosity (Benzi et al. 1998)

$$\frac{du_n}{dt}=i\left(k_{n+1}u_{n+1}^{}u_{n+2}\frac{1}{2}k_nu_{n1}^{}u_{n+1}+\frac{1}{2}k_{n1}u_{n2}u_{n1}\right)\nu _n^{(e)}k_n^2u_n+f_n$$
(22)   
where now $`n=1,\mathrm{},N_M<N`$ and the eddy viscosity, restricted to the last two shells, has the form

$$\nu _n^{(e)}=\kappa \frac{|u_n|}{k_n}\left(\delta _{n,N_M1}+\delta _{n,N_M}\right)$$
(23)   
where $`\kappa `$ is a constant of order $`1`$ (see Appendix C). The model equations (22) are the analogous of large eddy simulation (LES) in Shell Model which is one of the most popular numerical method for integrating large scale flows. Thus, although Shell Models are not realistic models for large scale geophysical flows (being nevertheless a good model for small scale turbulent fluctuations), the study of the effect of truncation in term of eddy viscosity is of general interest.

In Figure 5 we show $`\lambda _{MM}(\delta )`$, i.e. the FSLE computed for the model equations (22) with $`N=24`$ at different resolutions $`N_M=9,15,20`$. A plateau is detected for small amplitudes of the error $`\delta `$, corresponding to the leading Lyapunov exponent, which increases with increasing resolution – being proportional to the fastest timescale – according to $`\lambda k_{N_M}^{2/3}`$. At larger $`\delta `$ the curves collapse onto the $`\lambda _{TT}(\delta )`$, showing that large-scale statistics of the model is not affected by the small-scales resolution.

The capability of the model to predict satisfactorily the statistical features of the “true” dynamics is not anyway determined by $`\lambda _{MM}(\delta )`$ but by $`\lambda _{TM}(\delta )`$, which is shown in Figure 6.

Increasing the resolution $`N_M=9,15,20`$ towards the fully resolved case $`N=24`$ the model improves, in agreement with the expectation that $`\lambda _{TM}`$ approaches $`\lambda _{TT}`$ for a perfect model. At large $`\delta `$ the curves practically coincide, showing that the predictability time for large error sizes (associated with large scales) is independent on the details of small-scale modeling. Better resolved models achieve $`\lambda _{TM}\lambda _{TT}`$ for smaller values of the error $`\delta `$.

## V Conclusions

In this Paper the effects of the uncertainty of the evolution laws on the predictability properties are investigated and quantitatively characterized by means of the Finite Size Lyapunov Exponent. In particular, we have considered systems involving several characteristic scales and times. In these cases, it is rather natural to investigate what is the effect of small scale parameterization on large scale dynamics.

It has been shown that in systems where there is a negligible feedback on the small scales by the large ones, the dynamics of the former ones can be thoroughly discarded, without affecting the statistical features of large scales and the ability to forecast them. On the other side, when this feedback is present, the crude approximation of cutting the small scale variables off is no longer acceptable. In this case one has to model the action of fast modes (small scales) on slow modes (large scales) with some effective term, in order to recover a satisfactory forecasting of large scales. The renowned eddy-viscosity modelization is an instance of the general modeling scheme that has been here discussed.

## VI Acknowledgments

We thank L. Biferale for useful suggestions and discussions. This work was partially supported by INFM (Progetto Ricerca Avanzata TURBO) and by MURST (program 9702265437). A special acknowledgment goes to B. Marani for warm and continuous support.

## A An example of structural unstable system

In order to see that a non generic perturbation, although very “small”, can produce dramatic changes in the dynamics, let us discuss a simple example following (Berkooz 1994; Holmes et al. 1996). We consider the one-dimensional chaotic map $`x_{t+1}=f(x_t)`$ with $`f(x)=4x`$ mod $`1`$, and a perturbed version of it:

$$f_p(x)=\{\begin{array}{cc}8x\frac{9}{2}\hfill & x[\frac{5}{8},\frac{247}{384}]\hfill \\ & \\ \frac{1}{2}x+\frac{1}{3}\hfill & x[\frac{247}{384},\frac{265}{384}]\hfill \\ & \\ 8x\frac{29}{6}\hfill & x[\frac{265}{384},\frac{17}{24}]\hfill \\ & \\ 4x\text{mod}\mathrm{\hspace{0.33em}1}\hfill & \text{otherwise}.\hfill \end{array}$$
(A1)   
The perturbed map is identical to the original outside the interval $`[5/8,17/24]`$, and the perturbation is very small in $`L_2`$ norm. Nevertheless, the fixed point $`x=\frac{2}{3}`$, which is unstable in the original dynamics, becomes stable in the perturbed one. Moreover it is a global attractor for $`f_p(x)`$, i.e. almost every point in $`[0,1]`$ asymptotically approaches $`x=\frac{2}{3}`$ (see Figure 1).

Now, if one compares the trajectories obtained iterating $`f(x)`$ or $`f_p(x)`$ it is not difficult to understand that orbits starting outside $`[5/8,17/24]`$ remain identical for a certain time but unavoidably they differ utterly in the long time behavior. It is easy to realize that the transient chaotic behavior of the perturbed orbits can be rendered arbitrarily long by reducing the interval in which the two dynamics differ. This example shows how even an ostensibly small perturbation (in usual norms) can modify dramatically the dynamics.

## B Computation of the Finite size Lyapunov exponent

In this appendix we discuss in detail the computation of the Finite Size Lyapunov Exponent for both continuous dynamics (differential equations) and discrete dynamics (maps).

The practical method for computing the FSLE goes as follows. Defined a given norm for the distance $`\delta (t)`$ between the reference and perturbed trajectories, one has to define a series of thresholds $`\delta _n=r^n\delta _0`$ ($`n=1,\mathrm{},N`$), and to measure the “doubling times” $`T_r(\delta _n)`$ that a perturbation of size $`\delta _n`$ takes to grow up to $`\delta _{n+1}`$. The threshold rate $`r`$ should not be taken too large, because otherwise the error has to grow through different scales before reaching the next threshold. On the other hand, $`r`$ cannot be too close to one, because otherwise the doubling time would be of the order of the time step in the integration. In our examples we typically use $`r=2`$ or $`r=\sqrt{2}`$. For simplicity $`T_r`$ is called “doubling time” even if $`r2`$.

The doubling times $`T_r(\delta _n)`$ are obtained by following the evolution of the separation from its initial size $`\delta _{min}\delta _0`$ up to the largest threshold $`\delta _N`$. This is done by integrating the two trajectories of the system starting at an initial distance $`\delta _{min}`$. In general, one must choose $`\delta _{min}\delta _0`$, in order to allow the direction of the initial perturbation to align with the most unstable direction in the phase-space. Moreover, one must pay attention to keep $`\delta _N<\delta _{saturation}`$, so that all the thresholds can be attained ($`\delta _{saturation}`$ is the typical distance of two uncorrelated trajectory, i.e. the size of the attractor). For the second kind predictability problem, i.e. the computation of $`\lambda _{TM}(\delta )`$, one can safely take $`\delta _{min}=0`$ because this do not prevent the separation of trajectories.

The evolution of the error from the initial value $`\delta _{min}`$ to the largest threshold $`\delta _N`$ carries out a single error-doubling experiment. At this point one rescales the model trajectory at the initial distance $`\delta _{min}`$ with respect to the true trajectory and starts another experiment. After $`𝒩`$ error-doubling experiments, we can estimate the expectation value of some quantity $`A`$ as:

$$A_e=\frac{1}{𝒩}\underset{i=1}{\overset{𝒩}{}}A_i.$$
(B1)   
This is not the same as taking the time average as in (6) because different error doubling experiments may takes different times. Indeed we have

$$A_t=\frac{1}{T}_0^TA(t)𝑑t=\frac{\underset{i}{}A_i\tau _i}{_i\tau _i}=\frac{A\tau _e}{\tau _e}.$$
(B2)   
In the particular case in which $`A`$ is the doubling time itself we have from (6) and (B2)

$$\lambda (\delta _n)=\frac{1}{T_r(\delta _n)_e}\mathrm{ln}r.$$
(B3)   
The method described above assumes that the distance between the two trajectories is continuous in time. This is not true for maps of for discrete sampling in time and the method has to be slightly modified. In this case $`T_r(\delta _n)`$ is defined as the minimum time at which $`\delta (T_r)r\delta _n`$. Because now $`\delta (T_r)`$ is a fluctuating quantity, from (B2) we have

$$\lambda (\delta _n)=\frac{1}{T_r(\delta _n)_e}\mathrm{ln}\left(\frac{\delta (T_r)}{\delta _n}\right)_e,.$$
(B4)   
We conclude by observing that the computation of the FSLE is not more expensive than the computation of the Lyapunov exponent by standard algorithm. One has simply to integrate two copies of the system (or two different systems for second kind predictability) and this can be done also for very complex simulations.

## C Parameterization of small scales

Typically a realistic problem (e.g. turbulence) involves many interacting degrees of freedom with different characteristic times. Let us indicate with $`𝐳`$ the state of the system under consideration, with an evolution law:

$$\frac{\mathrm{d}𝐳}{\mathrm{d}t}=𝐅(𝐳),𝐅,𝐳^N.$$
(C1)   
The dynamical variables $`𝐳`$ can be split in two sets:

$$𝐳=(𝐱,𝐲),$$
(C2)   
where $`𝐱^n`$ and $`𝐲^m`$ ($`N=n+m`$), being respectively $`𝐱`$ and $`𝐲`$ the “slow” and “fast” variables. The distinction between slow and fast variables is often largely arbitrary.

The evolution equation (C1) is divided into two blocks, the first one containing the dynamics of the slow variables, the second one associated with the dynamics of the fast variables

$$\{\begin{array}{c}\frac{d𝐱}{dt}=𝐅_\mathrm{𝟏}(𝐱)+𝐅_\mathrm{𝟐}(𝐱,𝐲)\\ \frac{d𝐲}{dt}=\stackrel{~}{𝐅}_\mathrm{𝟏}(𝐱,𝐲)+\stackrel{~}{𝐅}_\mathrm{𝟐}(𝐲)\end{array}$$
(C3)   
If one is interested only in the slow variables it is necessary to write an “effective” equation for $`𝐱`$. As far as we know there is only one case for which it is simple to find the effective equations for $`𝐱`$. If the characteristic times of the fast variables are much smaller than those ones of the $`𝐱`$ (adiabatic limit), one can write:

$$𝐲=<𝐲>+𝜼(t)$$
(C4)   
where $`𝜼`$ is a Wiener process, i.e. a zero mean Gaussian process with

$$<\eta _i(t)\eta _j(t^{^{}})>=<\delta y_i^2>\delta _{ij}\delta (tt^{^{}}).$$
(C5)   
Therefore one obtains for the slow variables:

$$\frac{\mathrm{d}𝐱}{\mathrm{d}t}=𝐅_\mathrm{𝟏}(𝐱)+\delta 𝐅_\mathrm{𝟏}(𝐱)+\delta 𝐖(𝐱,𝜼)$$
(C6)   
where $`\delta 𝐅_\mathrm{𝟏}(𝐱)=𝐅_\mathrm{𝟐}(𝐱,<𝐲>)+\delta 𝐅_\mathrm{𝟐}`$, $`\delta F_{2,j}=1/2_i^2F_{2,j}/y_j^2<\delta y_i^2>`$ and $`\delta W_i=_iF_{2,j}/y_j|_{<𝐲>}\eta _i(t)`$. Basically the slow variables $`𝐱`$ obey to a non linear Langevin equation.

Here the role of the fast degrees of freedom becomes relatively simple: they give small changes to the drift $`𝐅_\mathrm{𝟏}𝐅_\mathrm{𝟏}+\delta 𝐅_\mathrm{𝟏}`$ and a noise term $`\delta 𝐖(𝐱,𝜼)`$. We remark that the validity of the above argument is rather limited. Even if one has a large time scale separation, the statistics of the fast variables can be very far from the Gaussian distribution. In particular, in system with feedback ($`\stackrel{~}{𝐅}_10`$) one cannot model the fast variable $`𝐲`$ independently of the resolved $`𝐱`$.

In the generic situation the construction of the effective equation for $`𝐱`$ requires to follow phenomenological arguments which depend on the physical mechanism of the particular problem. For example, for the Lorenz ’96 model discussed in sect. III, where $`F_{2,k}(𝐱,𝐲)=_{j=1,J}y_{j,k}`$, we use the following procedure for the parameterization of the fast variables and the building of the effective eq. for $`𝐱`$. Instead of assuming (C4) we mimic the fast variables in terms of the slow ones:

$$𝐲(t)=𝐠(𝐱(t))=<𝐲|𝐱(t)>+\eta (t)$$
(C7)   
where $`<|𝐱>`$ stands for the conditional average and $`\eta (t)`$ is a noise term. Inserting (C7) into the first of (C3) one obtains

$$\frac{\mathrm{d}𝐱}{\mathrm{d}t}=𝐅_1(𝐱)+𝐅_2(𝐱,𝐲)=𝐅_1(𝐱)+𝐅_2(𝐱,<𝐲|𝐱>)+\delta 𝐅_2(𝐱)$$
(C8)   
where

$$\delta F_{2,i}=\underset{j,k}{}\frac{^2F_{2,i}}{y_jy_k}|_{y=<y|x>}\eta _j\eta _k$$
(C9)   
In the Lorenz ’96 model (14), because of the linear coupling between the different scales, the terms $`\delta 𝐅_2`$ are absent and one has a close model for the large scale variables

$$\frac{\mathrm{d}𝐱}{\mathrm{d}t}=𝐅_\mathrm{𝟏}(𝐱)+𝐅_\mathrm{𝟐}(𝐱,<𝐲|𝐱>)$$
(C10)   
The ansatz (C7) is well verified in the numerical simulations. We have computed the $`\lambda _{TM}(\delta )`$ by using a best fit for $`𝐅_\mathrm{𝟐}`$ and we have obtained a good reproduction of the $`\lambda _{TT}(\delta )`$ for large $`\delta `$. In the Lorenz ’96 model (14), where the coupling between slow and fast variables is practically linear, one has that $`F_{2,k}(𝐱,<𝐲|𝐱>)=_{j=1,J}<y_{j,k}|x_k>\nu _ex_k`$.

Now we will discuss the case of the Shell Model parameterization which pertains to the general issue of the subgrid-scale modelization. The literature on this field and the related problems (e.g. closure in fully developed turbulence) is enormous and we do not pretend to discuss here in details this field. Let us only recall the basic idea introduced over a century ago by Boussinesq, and later developed further by Taylor, Prandtl and Heisenberg – to cite some of the most famous ones– for fully developed turbulence (Frisch 1995). In a nutshell the idea is to mimic the energy flux from the large to the small scales (in our terms from slow to fast variables) by an effective dissipation: the effect of the small scales on the large ones can be modeled as an enhanced molecular viscosity.

By simple dimensional arguments one can argue that the effects of small scales can be replaced by an effective viscosity at scales $`r`$, given by

$$\nu ^{(e)}r\delta v(r)$$
(C11)   
where $`\delta v(r)`$ is the velocity fluctuation on the scale $`r`$.

The above argument for the Shell Model (17) gives (Benzi et al. 1998):

$$\nu _n^{(e)}=\kappa \frac{|u_n|}{k_n}$$
(C12)   
where $`\kappa O(1)`$ is an empirical constant. From eq. (C11) one could naively think to use dimensional argument á la Kolmogorov to set a constant eddy viscosity $`\nu _n^{(e)}k_n^{4/3}`$. In this way one forgets the dynamics and this can cause numerical blow up. More sophisticated arguments that do not include the dynamics lead to similar problems.

Let us remark that the parameterization (C12) is not exactly identical to those obtained by closure approaches where the eddy viscosity is given in terms of averaged quantities. In our case this would mean to write $`u_n^2^{1/2}`$ instead of $`|u_n|`$ in (C12).

After this discussion it is easy to recognize that the parameterization in terms of conditional averages introduced for the Lorenz ’96 model is, a posteriori, an eddy viscosity model.

FIGURE CAPTIONS

* The map $`f_p`$ of equation (A1) (solid line) and the original chaotic map $`f`$ (dashed line).
* Finite Size Lyapunov Exponents $`\lambda _{TT}(\delta )`$ ($`+`$) and $`\lambda _{TM}(\delta )`$ ($`\times `$) versus $`\delta `$ for the Lorenz model (4) with $`\sigma =c=10`$, $`b=8/3`$, $`R=45`$ and $`ϵ=0.001`$. The dashed line represents the leading Lyapunov exponent for the unperturbed system ($`\lambda 1.2`$). The statistics is over $`10^4`$ realizations.
* $`\lambda _{TT}(\delta )`$ (crosses, $`\times `$) and $`\lambda _{TM}(\delta )`$ (open squares, $`\mathrm{}`$) versus $`\delta `$ for the Rayleigh-Bénard model (7) with $`C=0.5`$, $`k_L=1`$, $`\omega _L=1`$, $`B_L=0.3`$, $`k_S=4`$, $`\omega _S=4`$, $`B_S=0.3`$ and $`ϵ=0.125`$. The straight line indicates the $`\delta ^2`$ slope. The statistics is over $`10^4`$ realizations.
* Finite Size Lyapunov Exponents for the Lorenz ’96 model $`\lambda _{TT}(\delta )`$ (solid line) and $`\lambda _{TM}(\delta )`$ versus $`\delta `$ obtained by dropping the fast modes ($`+`$) and with eddy viscosity parameterization ($`\times `$) as discussed in (15) and (16). The parameters are $`F=10`$, $`K=36,J=10`$, $`\nu =1`$ and $`c=b=10`$, implying that the typical $`y`$ variable is $`10`$ times faster and smaller than the $`x`$ variable. The value of the parameter $`\nu _e=4`$ is chosen after a numerical integration of the complete equations as discussed in Appendix C. The statistics is over $`10^4`$ realizations.
* The FSLE for the eddy-viscosity shell model (22) $`\lambda _{MM}(\delta )`$ at various resolutions $`N_M=9(+),15(\times ),20()`$. For comparison it is drawn the FSLE $`\lambda _{TT}(\delta )`$ (continuous line). Here $`\kappa =0.4`$, $`k_0=0.05`$.
* The FSLE between the eddy-viscosity shell model and the full shell model $`\lambda _{TM}(\delta )`$, at various resolutions $`N_M=9(+),15(\times ),20()`$. For comparison it is drawn the FSLE $`\lambda _{TT}(\delta )`$ (continuous line). The total number of shell for the complete model is $`N=24`$, with $`k_0=0.05`$, $`\nu =10^7`$.