# Non Asymptotic Properties of Transport and Mixing

## I Introduction

Transport processes play a crucial role in many natural phenomena. Among the many examples, we just mention the particle transport in geophysical flows which is of obvious interest for atmospheric and oceanic issues. The most natural framework for investigating such phenomena is to adopt a Lagrangian viewpoint in which the particles are advected by a given Eulerian velocity field $`𝒖(𝒙,t)`$ according to the differential equation

$$\frac{d𝒙}{dt}=𝒖(𝒙,t)=𝒗(t),$$
(1)   
where, by definition, $`𝒗(t)`$ is the Lagrangian particle velocity.

Despite the apparent simplicity of (1), the problem of connecting the Eulerian properties of $`𝒖`$ to the Lagrangian properties of the trajectories $`𝒙(t)`$ is a very difficult task. In the last $`2030`$ years the scenario has become even more complex by the recognition of the ubiquity of Lagrangian chaos (chaotic advection). Even very simple Eulerian fields can generate very complex Lagrangian trajectories which are practically indistinguishable from those obtained in a complex, turbulent, flow .

Despite these difficulties, the study of the relative dispersion of two particles can give some insight on the link between Eulerian and Lagrangian properties at different length-scales. Indeed, the evolution of the separation $`𝑹(t)=𝒙^{(2)}(t)𝒙^{(1)}(t)`$ between two tracers is given by

$$\frac{d𝑹}{dt}=𝒗^{(2)}(t)𝒗^{(1)}(t)=𝒖(𝒙^{(1)}(t)+𝑹(t),t)𝒖(𝒙^{(1)}(t),t)$$
(2)   
and thus depends on the velocity difference on scale $`𝑹`$. It is obvious from (2) that Eulerian velocity components of typical scale much larger than $`𝑹`$ will not contribute to the evolution of $`𝑹`$. Since, in incompressible flows, separation $`𝑹`$ typically grows in time  we have the nice situation in which from the evolution of the relative separation we can, in principle, extract the contributions of all the components of the velocity field. For these reasons, in this paper we prefer to study relative dispersion instead of absolute dispersion. For spatially infinite cases, without mean drift there is no difference; for closed basins the relative dispersion is, for many aspects, more interesting than the absolute one, which is dominated by the sweeping induced by large scale flow.

There are very few general results on the link between Eulerian and Lagrangian properties and only for asymptotic behaviors. Let us suppose that the Eulerian velocity field is characterized by two typical length-scales: the (small) scale $`l_u`$ below which the velocity is smooth, and a (large) scale $`L_0`$ representing the size of the largest structures present in the flow. Of course, in most non turbulent flows it will turn out that $`ł_uL_0`$.

At very small separations $`Rł_u`$ we have that the velocity difference in (2) can be reasonably approximated by a linear expansion in $`R`$, which in most time-dependent flows leads to an exponential growth of the separation of initially close particles, a phenomenon known as Lagrangian chaos

$$\mathrm{ln}R(t)\mathrm{ln}R(0)+\lambda t$$
(3)   
(the average is taken over many couples with initial separation $`R(0)`$). The coefficient $`\lambda `$ is the Lagrangian Lyapunov exponent of the system . The rigorous definition of the Lyapunov exponent imposes to take the two limits $`R(0)0`$ and then $`t\mathrm{}`$: in physical terms these limits amount to the requirement that the separation has not to exceed the scale $`ł_u`$ but for very large times. This is a very strict condition, rarely accomplished in real flows, rendering often infeasible the experimental observation of the behavior (3).

On the opposite limit, for very long times and for separations $`RL_0`$, the two trajectories $`𝒙^{(1)}(t)`$ and $`𝒙^{(2)}(t)`$ feel two velocities which can practically be considered as uncorrelated. We thus expect normal diffusion, i.e.

$$R^2(t)2Dt.$$
(4)   
Also in this case it is necessary to remark that the asymptotic behavior (4) cannot be attained in many realistic situations, the most common of which is the presence of boundaries at a scale comparable with $`L_0`$. In absence of boundaries it is possible to formulate sufficient conditions on the nature of the Eulerian flow, under which normal diffusion (4) always takes place asymptotically .

Between the two asymptotic regimes (3) and (4) the behavior of $`R(t)`$ depends on the particular flow. The study of the evolution of the relative dispersion in this crossover regime is very interesting and can give an insight on the Eulerian structure of the velocity field.

To summarize, in all systems in which the characteristic length-scales are not sharply separated, it is not possible to describe dispersion in terms of asymptotic quantities. In such cases, different approaches are required. Let us mention some examples: the symbolic dynamics approach to the sub-diffusive behavior in a stochastic layer  and to mixing in meandering jets ; the study of tracer dynamics in open flows in terms of chaotic scattering  and the exit time description for transport in semi-enclosed basins  and open flows .

The aim of the present paper is to discuss the use of an indicator – the Finite Size Lyapunov Exponents (FSLE), originally introduced in the context of predictability problems  – to study and characterize non-asymptotic transport in non-ideal systems, e.g. closed basins and systems in which the characteristic length-scales are not sharply separated.

In section II we introduce the basic tools for the finite-scale analysis and we discuss their general properties. Section III is devoted to the evaluation of our method on some numerical examples. We shall see that even in apparently simple situations the use of finite scale analysis avoids possible misinterpretation of the results. In section IV the method is applied to two physical problems: the analysis of experimental drifter data and the numerical study of relative dispersion in fully developed turbulence. Conclusions are presented in section V. The appendices report, for sake of self-consistency, some technical aspects.

## II Finite size diffusion coefficient

In order to introduce the finite size analysis for the dispersion problem let us start with a simple example. We consider a set of N particle pairs advected by a smooth (e.g. spatially periodic) velocity field with characteristic length $`l_u`$. Denoting with $`R_i^2(t)`$ the square separation of the $`i`$-th couple, we define

$$R^2(t)=\frac{1}{N}\underset{i=0}{\overset{N}{}}R_i^2.$$
(5)   
We assume that the Lagrangian motion is chaotic, thus we expect the following regimes to hold

$$R^2(t)\{\begin{array}{cc}R_0^2\mathrm{exp}(L(2)t)\hfill & \text{if }R^2(t)^{1/2}l_u\hfill \\ 2Dt\hfill & \text{if }R^2(t)^{1/2}l_u\hfill \end{array},$$
(6)   
where $`L(2)2\lambda `$ is the generalized Lyapunov exponent , $`D`$ is the diffusion coefficient and we assume that $`R_i(0)=R_0`$.

An alternative method to characterize the dispersion properties is by introducing the “doubling time” $`\tau (\delta )`$ at scale $`\delta `$ as follows : given a series of thresholds $`\delta ^{(n)}=r^n\delta ^{(0)}`$, one can measure the time $`T_i(\delta ^{(0)})`$ it takes for the separation $`R_i(t)`$ to grow from $`\delta ^{(0)}`$ to $`\delta ^{(1)}=r\delta ^{(0)}`$, and so on for $`T_i(\delta ^{(1)}),T_i(\delta ^{(2)}),\mathrm{}`$ up to the largest considered scale. The $`r`$ factor may be any value $`>\mathrm{\hspace{0.17em}1}`$, properly chosen in order to have a good separation between the scales of motion, i.e. $`r`$ should be not too large. Strictly speaking, $`\tau (\delta )`$ is exactly the doubling time if $`r=\mathrm{\hspace{0.17em}2}`$.

Performing the doubling time experiments over the $`N`$ particle pairs, one defines the average doubling time $`\tau (\delta )`$ at scale $`\delta `$ as

$$\tau (\delta )=<T(\delta )>_e=\frac{1}{N}\underset{i=1}{\overset{N}{}}T_i(\delta ).$$
(7)   
It is worth to note that the average (7) is different from the usual time average (see Appendix A).

Now we can define the Finite Size Lagrangian Lyapunov Exponent (see  for a detailed discussion) in terms of the average doubling time as

$$\lambda (\delta )=\frac{\mathrm{ln}r}{\tau (\delta )},$$
(8)   
which quantifies the average rate of separation between two particles at a distance $`\delta `$. Let us remark that $`\lambda (\delta )`$ is independent of $`r`$, for $`r`$ close to $`1`$. For very small separations (i.e. $`\delta l_u`$) one recovers the Lagrangian Lyapunov exponent $`\lambda `$,

$$\lambda =\underset{\delta 0}{lim}\frac{1}{\tau (\delta )}\mathrm{ln}r.$$
(9)   
In this framework the finite size diffusion coefficient , $`D(\delta )`$, dimensionally turns out to be

$$D(\delta )=\delta ^2\lambda (\delta ).$$
(10)   
Note the absence of the factor $`2`$, as one may expect from (6), in the denominator of $`D(\delta )`$; this is because $`\tau (\delta )`$ is a difference of times. For a standard diffusion process $`D(\delta )`$ approaches the diffusion coefficient $`D`$ (see eq. (6)) in the limit of very large separations ($`\delta l_u`$). This result stems from the scaling of the doubling times $`\tau (\delta )\delta ^2`$ for normal diffusion.

Thus the finite size Lagrangian Lyapunov exponent $`\lambda (\delta )`$ behaves as follows:

$$\lambda (\delta )\{\begin{array}{cc}\lambda \hfill & \text{if }\delta l_u\hfill \\ D/\delta ^2\hfill & \text{if }\delta l_u\hfill \end{array},$$
(11)   
One could naively conclude, matching the behaviors at $`\delta l_u`$, that $`D\lambda l_u^2`$. This is not always true, since one can have a rather large range for the crossover due to nontrivial correlations which can be present in the Lagrangian dynamics .

One might wonder that the introduction of $`\tau (\delta )`$ is just another way to look at $`R^2(t)`$. This is true only in limiting cases, when the different characteristic lengths are well separated and intermittency is weak. A similar idea of using times for the computation of the factor diffusion coefficient in nontrivial cases was developed in Ref. .

If one wants to identify the physical mechanisms acting on a given spatial scale, the use of scale dependent quantities is more appropriate than time dependent ones.

For instance, in presence of strong intermittency (which is indeed a rather usual situation) $`R^2(t)`$ as a function of $`t`$ can be very different in each realization. Typically one has (see figure 1a), different exponential growth rates for different realizations, producing a rather odd behavior of the average $`R^2(t)`$ not due to any physical mechanisms. For instance in figure 1b we show the average $`R^2(t)`$ versus time $`t`$; at large times one recovers the diffusive behavior but at intermediate times appears an “anomalous” diffusive regime which is only due to the superposition of exponential and diffusive contributions by different samples at the same time. On the other hand, by exploiting the tool of doubling times one has an unambiguous result (see figure 1c) .

An important physical problem where the behavior of $`\tau (\delta )`$ is essentially well understood is the relative dispersion in 3D fully developed turbulence. Here the smallest Eulerian scale $`l_u`$ is the Kolmogorov scale at which the flow becomes smooth. In the inertial range $`l_u<R<L_0`$ we expect the Richardson law to hold $`R^2(t)t^3`$; for separations larger than the integral scale $`L_0`$ we have normal diffusion. In terms of the finite size Lyapunov exponent we thus expect three different regimes:

1. $`\lambda (\delta )=\lambda `$ for $`\delta l_u`$
2. $`\lambda (\delta )\delta ^{2/3}`$ for $`l_u\delta L_0`$
3. $`\lambda (\delta )\delta ^2`$ for $`\delta L_0`$

We will see in section IV than even for large Reynolds numbers, the characteristic lengths $`l_u`$ and $`L_0`$ are not sufficiently separated and the different scaling regimes for $`R^2(t)`$ cannot be well detected. The fixed scale analysis in terms of $`\lambda (\delta )`$ for fully developed turbulence presents clear advantages with respect to the fixed time approach.

## III Numerical Results on simple flows

In this section we shall discuss some examples of applications of the above introduced indicator $`\lambda (\delta )`$ (or equivalently $`D(\delta )`$) for simple flows. The technical and numerical details of the finite size Lyapunov exponent computation are settled out in Appendix A.

In a generic case in addition to the two asymptotic regimes (11) discussed in section II, we expect another universal regime due to the presence of the boundary of given size $`L_B`$. For separations close to the saturation value $`\delta _{max}L_B`$ we expect the following behavior to hold for a broad class of systems :

$$\lambda (\delta )=\frac{D(\delta )}{\delta ^2}\frac{(\delta _{max}\delta )}{\delta }.$$
(12)   
The proportionality constant is given by the second eigenvalue of the Perron-Frobenius operator which is related to the typical time of exponential relaxation of tracers’ density to uniform distribution (see Appendix B).

### A A model for transport in Rayleigh-Bénard convection

The advection in two dimensional incompressible flows in absence of molecular diffusion is given by Hamiltonian equation of motion where the stream function, $`\psi `$, plays the role of the Hamiltonian:

$$\frac{dx}{dt}=\frac{\psi }{y},\frac{dy}{dt}=\frac{\psi }{x}.$$
(13)   
If $`\psi `$ is time-dependent one typically has chaotic advection. As an example let us consider the time-periodic Rayleigh-Bénard convection, which can be described by the following stream function :

$$\psi (x,y,t)=\frac{A}{k}\mathrm{sin}\left\{k\left[x+B\mathrm{sin}(\omega t)\right]\right\}W(y),$$
(14)   
where $`W(y)`$ satisfies rigid boundary conditions on the surfaces $`y=0`$ and $`y=a`$ (we use $`W(y)=\mathrm{sin}(\pi y/a)`$). The two surfaces $`y=a`$ and $`y=0`$ are the top and bottom surfaces of the convection cell. The time dependent term $`B\mathrm{sin}(\omega t)`$ represents lateral oscillations of the roll pattern which mimic the even oscillatory instability .

Concerning the analysis in terms of the finite size Lyapunov exponent one has that, if $`\delta `$ is much smaller than the domain size, $`\lambda (\delta )=\lambda `$. At larger values of $`\delta `$ we find standard diffusion $`\lambda (\delta )=D/\delta ^2`$ with good quantitative agreement with the value of the diffusion coefficient evaluated by the standard technique, i.e. using $`R^2(t)`$ as a function of time $`t`$.

In order to study the effects of finite boundaries on the diffusion properties we confine the tracers’ motion in a closed domain. This can be achieved by slightly modifying the stream function (14). We have modulated the oscillating term in such a way that for $`|x|=L_B`$ the amplitude of the oscillation is zero, i.e. $`BB\mathrm{sin}(\pi x/L_B)`$ with $`L_B=2\pi n/k`$ ($`n`$ is the number of convective cells). In this way the motion is confined in $`x[L_B,L_B]`$.

In figure 2 we show $`\lambda (\delta )`$ for two values of $`L_B`$. If $`L_B`$ is large enough one can distinguish the three regimes: exponential, diffusive and the saturation regime eq. (12). Decreasing the size of the boundary $`L_B`$, the range of the diffusive regime decreases, while for small values of $`L_B`$, it disappears.

### B Point vortices in a Disk

We now consider a two-dimensional time-dependent flow generated by $`M`$ point vortices, with circulations $`\mathrm{\Gamma }_1,\mathrm{},\mathrm{\Gamma }_M`$, in a disk of unit radius . The passive tracers are advected by the time dependent velocity field generated by the vortices and behave chaotically for any $`M>2`$. Let us note that in this case the scale separation is not imposed by hand, but depends on $`M`$ and on the energy of the vortex system . Figure 3a shows the relative diffusion as a function of time in a system with $`M=4`$ vortices. Apparently there is an intermediate regime of anomalous diffusion. However from figure 3b one can clearly see that, with the fixed scale analysis, only two regimes survive: exponential and saturation. Comparing figure 3a and figure 3b one understands that the appearance of the spurious anomalous diffusion regime in the fixed time analysis is due to the mechanism described in section II.

The absence of the diffusive regime $`\lambda (\delta )\delta ^2`$ is due to the fact that the characteristic length of the velocity field, which is comparable with the typical distance between two close vortices, is not much smaller than the size of the basin.

### C Random walk on a fractal object: an anomalous diffusive case

In this section we discuss the case of particles performing a continuous random walk on a fractal object of fractal dimension $`D_F`$, where one has sub-diffusion. We show that also in situation of anomalous diffusion (e.g. sub-diffusion) the FSLE is able to recognize the correct behavior. In the following section we consider the case of fully developed turbulence which displays super-diffusion.

In a fractal object due to the presence of voids, i.e. forbidden regions for the particles, one expects a decreasing of the spreading, and because of the self-similar structure of the domain (i.e. voids on all scales) a sub-diffusive behavior is expected. It is worth to note that the particles do not diffuse with the same law from any points of the domain (due to the presence of voids), hence in order to define a diffusive-like behavior one has to average over all possible particles’ position. For discrete random walk on a fractal lattice it is known that the diffusion follows the law $`<R^2(t)>t^{2/D_W}`$ with $`D_W>2`$, i.e. sub-diffusion . The quantity $`D_W`$ is related to the spectral or fracton dimension $`D_S`$, by the relation $`D_W=2D_F/D_S`$, and it depends on the detailed structure of the fractal object .

We study the relative dispersion of 2-D continuous random walk in a Sierpinsky Carpet with fractal dimension $`D_F=\mathrm{log}\mathrm{\hspace{0.17em}8}/\mathrm{log}\mathrm{\hspace{0.17em}3}`$. In our computation we use a resolution $`3^5`$, i.e. the fractal is approximated by five steps of the recursive building rule, in practice we perform a continuous random walk in a basin obtained with the above approximation of the Sierpinsky Carpet. We initialize the particles inside one of the smallest resolved structures, then we follow the growth of the relative dispersion with the FSLE method, and redeploy the particles in a small cell randomly chosen at the beginning of each doubling time experiment. From fig. 4 one can see that $`\lambda (\delta )\delta ^{1/.45}`$ which is an indication of sub-diffusion the exponent is in good agreement with the usual relative dispersion analysis (see the inset of fig. 4).

## IV Application of the FSLE

### A Drifter in the Adriatic Sea: data analysis and modelization

Lagrangian data recorded within oceanographic programs in the Mediterranean Sea  offer the opportunity to apply the fixed scale analysis to a geophysical problem, for which the standard characterization of the dispersion properties gives poor information.

The Adriatic Sea is a semi-enclosed basin, about 800 by 200 $`km`$ wide, connected to the whole Mediterranean Sea through the Otranto Strait . We adopt the reference frame in which the $`x,y`$axes are aligned, respectively, with the short side (transverse direction), orthogonal to the coasts, and the long side (longitudinal direction), along the coasts.

We have computed relative dispersion along the two axes, $`R_x^2(t)`$, $`R_y^2(t)`$ and FSLE $`\lambda (\delta )`$. The number of selected drifters for the analysis is 37, distributed in 5 different deployments in the Strait of Otranto, happened during the period December 1994 - March 1997, containing respectively 4, 9, 7, 7 and 10 drifters. To get as high statistics as possible, even to the cost of losing information on the seasonal variability, we shift the time tracks of all of the 37 drifters to $`tt_0`$, where $`t_0`$ is the time of the deployment, so that the drifters can be treated as a whole cluster. Moreover, to restrict the analysis only to the Adriatic basin, we discarded a drifter as soon as its latitude goes south of $`39.5`$ N or its longitude goes beyond $`19.5`$ E.

Before presenting the results of the data analysis, let us introduce a simplified model for the Lagrangian tracers motion in the Adriatic Sea. We assume as main features of the surface circulation the following elements : the drifter motion is basically two-dimensional; the domain is a quasi-closed basin; an anti-clockwise coastal current; two large cyclonic gyres; some natural irregularities in the Lagrangian motion induced by small scale structures. On the basis of these considerations, we introduce a deterministic chaotic model with mixing properties for the Lagrangian drifters. The stream function is given by the sum of three terms:

$$\mathrm{\Psi }(x,y,t)=\mathrm{\Psi }_0(x,y)+\mathrm{\Psi }_1(x,y,t)+\mathrm{\Psi }_2(x,y,t)$$
(15)   
defined as follows:

$$\mathrm{\Psi }_0(x,y)=\frac{C_0}{k_0}[\mathrm{sin}(k_0(y+\pi ))+\mathrm{cos}(k_0(x+2\pi ))]$$
(16)   
$$\mathrm{\Psi }_i(x,y,t)=\frac{C_i}{k_i}\mathrm{sin}(k_i(x+ϵ_i\mathrm{sin}(\omega _it)))\mathrm{sin}(k_i(y+ϵ_i\mathrm{sin}(\omega _it+\varphi _i))),(i=1,2),$$
(17)   
where $`k_i=\mathrm{\hspace{0.17em}2}\pi /\lambda _i`$, for $`i=0,1,2`$, $`\lambda _i`$’s are the wavelengths of the spatial structures of the flow; analogously $`\omega _j=\mathrm{\hspace{0.17em}2}\pi /T_j`$, for $`j=1,2`$, and $`T_j`$’s are the periods of the perturbations. In the non-dimensional expression of the equations, the length and time units have been set to $`200km`$ and $`7.5days`$, respectively.

The stationary term $`\mathrm{\Psi }_0`$ defines the boundary large scale circulation with positive vorticity. The contribution of $`\mathrm{\Psi }_1`$ contains the two cyclonic gyres and it is explicitly time-dependent through a periodic perturbation. The term $`\mathrm{\Psi }_2`$ gives the motion over scales smaller than the size of the large gyres and it is time-dependent as well. The zero-value isoline is defined as the boundary of the basin.

According to observation, we have chosen the parameters so that the velocity range is around $`\mathrm{\hspace{0.33em}0.3}ms^1`$; the length scales of the Eulerian structures are $`L_B\mathrm{\hspace{0.33em}1000}km`$ (coastal current), $`L_0\mathrm{\hspace{0.33em}200}km`$ (gyres) and $`l_u\mathrm{\hspace{0.33em}50}km`$ (vortices); the typical recirculation times, for gyres and vortices, are $``$ 1 month and $``$ 1 week, respectively; the oscillation periods are $`10days`$ (gyres) and $`2days`$ (vortices).

Let us discuss now the comparison between data and model results. The relative dispersion along the two directions of the basin, for data and model trajectories, are shown in figures 5a,b. The results for the model are obtained from the spreading of a cluster of $`10^4`$ initial conditions. When a particle reaches the boundary ($`\mathrm{\Psi }=0`$) it is eliminated.

For the diffusion properties, one cannot expect a scaling for $`R_{x,y}^2(t)`$ before the saturation regime, since the Eulerian characteristic lengths are not too small compared with the basin size. Indeed, we do not observe a power law behavior neither for the experimental data nor for the numerical model.

Let us stress that by opportunely fitting the parameters, we could obtain the model curves even closer to the experimental ones, but this would not be very meaningful since there is no clear theoretical expectation in a transient regime.

Let us now discuss the finite size Lyapunov exponent. The analysis of the experimental data has been averaged over the total number of couples out of 37 trajectories, under the condition that the evolution of the distance between two drifters is no longer followed when any of the two exits the Adriatic basin.

In fig. 6 we show the FSLE for data and model. In our case, as discussed above, we are far from asymptotical conditions, therefore we do not observe the scaling $`\lambda (\delta )\delta ^2`$.

The $`\lambda _M(\delta )`$ obtained from the minimal chaotic model (15-17) shows the typical step-like shape of a system with two characteristic time scales, and offers a scenario about how the FSLE of real trajectories may come out.

The relevant fact is that the large-scale Lagrangian features are well reproduced, at least at a qualitative level, by a relatively simple model. We believe that this agreement is not due to a particular choice of the model parameters, but rather to the fact that transport is mainly dominated by large scales whereas small scale details play a marginal role.

It is evident the major advantages of FSLE with respect to the usual fixed time statistics of relative dispersion: from the relative dispersion analysis of fig. 5 we are unable to recognize the underlying Eulerian structures, while the FSLE of fig. 6 suggests the presence of structures on different scales and with different characteristic times. In conclusion, the fixed scale analysis gives information for discriminating among different models for the Adriatic Sea.

### B Relative dispersion in fully developed turbulence

We consider now the relative dispersion of particles pairs advected by an incompressible, homogeneous, isotropic, fully developed turbulent field. The Eulerian statistics of velocity differences is characterized by the Kolmogorov scaling $`\delta v(r)r^{1/3}`$, in an interval of scales $`ł_urL_0`$, called the inertial range, $`l_u`$ is now the Kolmogorov scale. Due to the incompressibility of the velocity field particles will typically diffuse away from each other . For pair separations less than $`l_u`$ we have exponential growth of the separation of trajectories, typical of smooth flows, whereas at separations larger than $`L_0`$ normal diffusion takes place. In the inertial range the average pair separation is not affected neither by large scale components of the flow, which simply sweep the pair, nor by small scale ones, whose intensity is low and which act incoherently. Accordingly, the separation $`R(t)`$ feels mainly the action of velocity differences $`\delta v(R(t))`$ at scale $`R`$. As a consequence of the Kolmogorov scaling the separation grows with the Richardson law 

$$R^2(t)t^3.$$
(18)   
Non-asymptotic behavior takes place in such systems whenever $`l_u`$ is not much smaller than $`L_0`$, that is when the Reynolds number is not high enough. As a matter of facts even at very high Reynolds numbers, the inertial range is still insufficient to observe the scaling (18) without any ambiguity. On the other hand, we shall show that FSLE statistics is effective already a relatively small Reynolds numbers.

In order to investigate the problem of relative dispersion at various scale separations a practical tool is the use of synthetic turbulent fields. In fact, by means of stochastic processes it is possible to build a velocity field which reproduces the statistical properties of velocity differences observed in fully developed turbulence . In order to avoid the difficulties related to the presence of sweeping in the velocity field, we limit ourselves to a correct representation of two-point velocity differences. In this case, if one adopts the reference frame in which one of the two tracers is at rest at the origin (the so called a Quasi-Lagrangian frame of reference), the motion of the second particle is ruled by the velocity difference in this frame of reference, which has the the same single time statistics of the Eulerian velocity differences . The detailed construction of the synthetic Quasi-Lagrangian velocity field is presented in Appendix C.

In figure 7 we show the results of simulations of pair dispersion by the synthetic turbulent field with Kolmogorov scaling of velocity differences at Reynolds number $`Re10^6`$ . The expected super-diffusive regime (18) can be well observed only for huge Reynolds numbers (see also Ref. . To explain the depletion of scaling range for the relative dispersion let us consider a series of pair dispersion experiments, in which a couple of particles is released at a separation $`R_0`$ at time $`t=0`$. At a fixed time $`t`$, as customarily is done, we perform an average over all different experiments to compute $`R^2(t)`$. But, unless $`t`$ is large enough that all particle pairs have “forgotten” their initial conditions, the average will be biased. This is at the origin of the flattening of $`R^2(t)`$ for small times, which we can call a crossover from initial condition to self similarity. In an analogous fashion there is a crossover for large times, of the order of the integral time-scale, since some couples might have reached a separation larger than the integral scale, and thus diffuse normally, meanwhile other pairs still lie within the inertial range, biasing the average and, again, flattening the curve $`R^2(t)`$. This correction to a pure power law is far from being negligible for instance in experimental data where the inertial range is generally limited due to the Reynolds number and the experimental apparatus. For example, references  show quite clearly the difficulties that may arise in numerical simulations with the standard approach.

To overcome these difficulties we exploit the approach based on the fixed scale statistics. The outstanding advantage of averaging at a fixed separation scale is that it removes all crossover effects, since all sampled pairs belong to the inertial range. The expected scaling properties of the doubling times is obtained by a simple dimensional argument. The time it takes for particle separation to grow from $`R`$ to $`2R`$ can be estimate as $`T(R)R/\delta v(R)`$; we thus expect for the inverse doubling times the scaling

$$\frac{1}{T(R)}\frac{\delta v(R)}{R}R^{2/3}$$
(19)   
In figure 8 the great enhancement of the scaling range achieved by using the doubling times is evident. In addition, by using the FSLE it is possible to study in details the effect of Eulerian intermittency on the Lagrangian statistics of relative dispersion. See Ref. for a detailed discussion and a comparison with a multifractal scenario. The conclusion that can be drawn is that in this case doubling time statistics makes it possible a much better estimate of the scaling exponent with respect to the standard – fixed time – statistics.

## V Conclusions

In the study of relative dispersion of Lagrangian tracers one has to tackle situations in which the asymptotic behavior is never attained. This may happen in presence of many characteristic Eulerian scales or, what is typical of real systems, in presence of boundaries. It is worth to stress that such kind of systems are very common in geophysical flows , and also in plasma physics . Therefore a close understanding of non-asymptotic transport properties can give much relevant information about these natural phenomena.

To face these problems, in recent years, there have been proposed different approaches whose common ingredient is basically an “exit time” analysis. We remind the symbolic dynamics  and the chaotic scattering  approaches, the exit time description for transport in semi-enclosed basins , symplectic maps , open flows  and in plasma physics .

In this paper we have discussed the applications of the Finite Size Lyapunov Exponent, $`\lambda (\delta )`$, in the analysis of several situations. This method is based on the identification of the typical time $`\tau (\delta )`$ characterizing the diffusive process at scale $`\delta `$ through the exit time. This approach is complementary to the traditional one, in which one looks at the average size of the clouds of tracers as function of time. For values of $`\delta `$ much smaller than the smallest characteristic length of the Eulerian velocity field, one has that $`\lambda (\delta )`$ coincides with the maximum Lagrangian Lyapunov Exponent. For larger $`\delta `$ the shape of $`\lambda (\delta )`$ depends on the detailed mechanisms of spreading, i.e. the structure of the advecting velocity field and/or the presence of boundaries. The diffusive regime corresponds to the behavior $`\lambda (\delta )D/\delta ^2`$. If $`\delta `$ gets close to its saturation value, i.e. the characteristic size of the basin, the universal shape of $`\lambda (\delta )`$ can be obtained on the basis of dynamical system theory. In addition, we have shown that the fixed scale method is able to recognize the presence of a genuine anomalous diffusion.

A remarkable advantage of working at fixed scale (instead of at fixed time as in the traditional approach) is its ability to avoid misleading results, for instance apparent anomalous scaling over a certain time interval. Moreover, with the FSLE one obtains the proper scaling laws also for a relatively small inertial range for which the standard technique gives rather controversial answers.

The proposed method can be also applied in the analysis of drifter experimental data or in numerical model for Lagrangian transport.

## VI Acknowledgments

We thank V. Artale, E. Aurell, L. Biferale, P. Castiglione, A. Crisanti, M. Falcioni, R. Pasmanter, P.M. Poulain, M. Vergassola and E. Zambianchi for collaborations and discussions in last years. A particular acknowledgment to B. Marani for the continuous and warm encouragement. We are grateful to the ESF-TAO (Transport Processes in the Atmosphere and the Oceans) Scientific Program for providing meeting opportunities. This paper has been partly supported by INFM (Progetto di Ricerca Avanzato PRA-TURBO), MURST (no. 9702265437), and the European Network Intermittency in Turbulent Systems (contract number FMRX-CT98-0175).

## A Computation of the Finite size Lyapunov exponent

In this appendix we discuss in detail the method for computing 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{},P`$), 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 _P`$. 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 _P<\delta _{max}`$, so that all the thresholds can be attained ($`\delta _{max}`$ is the typical distance of two uncorrelated trajectory).

The evolution of the error from the initial value $`\delta _{min}`$ to the largest threshold $`\delta _P`$ 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 $`N`$ error-doubling experiments, we can estimate the expectation value of some quantity $`A`$ as:

$$A_e=\frac{1}{N}\underset{i=1}{\overset{N}{}}A_i.$$
(A1)   
This is not the same as taking the time average 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}.$$
(A2)   
In the particular case in which $`A`$ is the doubling time itself we have from (A2)

$$\lambda (\delta _n)=\frac{1}{T_r(\delta _n)_e}\mathrm{ln}r.$$
(A3)   
The method above described assumes that the distance between the two trajectories is continuous in time. This is not true for maps or for discrete sampling in time, thus 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 (A2) we have

$$\lambda (\delta _n)=\frac{1}{T_r(\delta _n)_e}\mathrm{ln}\left(\frac{\delta (T_r)}{\delta _n}\right)_e.$$
(A4)   
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 and this can be done also for very complex simulations.

## B Universal saturation behavior of $`\lambda (\delta )`$

In this appendix we present the derivation of the asymptotic behavior (12) of $`\lambda (\delta )`$ for $`\delta `$ close to the saturation. The computation is explicitly done for the simple case of a one dimensional Brownian motion in the domain $`[L_B,L_B]`$, with reflecting boundary conditions: the numerical simulations indicate that the result is of general applicability.

The evolution of the probability density $`p`$ is ruled by the Fokker-Planck equation

$$\frac{p}{t}=\frac{1}{2}D\frac{^2p}{x^2}$$
(B1)   
with the Neumann boundary conditions $`\frac{p}{x}(\pm L_B)=0`$ .

The general solution of (B1) is

$$p(x,t)=\underset{k=\mathrm{}}{\overset{\mathrm{}}{}}\widehat{p}(k,0)e^{ikx}e^{t/\tau _k}+c.c$$
(B2)   
where

$$\tau _k=\left(\frac{D}{2}\frac{\pi ^2}{L_B^2}k^2\right)^1,k=0,\pm 1,\pm 2,\mathrm{}$$
(B3)   
At large times $`p`$ approaches the uniform solution $`p_0=1/2L_B`$. Writing $`p`$ as $`p(x,t)=p_0+\delta p(x,t)`$ we have, for $`t\tau _1`$ ,

$$\delta p\mathrm{exp}(t/\tau _1).$$
(B4)   
The asymptotic behavior for the relative dispersion $`R^2(t)`$ is

$$R^2(t)=\frac{1}{2}(xx^{})^2p(x,t)p(x^{},t)𝑑x𝑑x^{}$$
(B5)   
For $`t\tau _1`$ using (B4) we obtain $`R^2(t)\left(\frac{L_B^2}{3}Ae^{t/\tau _1}\right)`$. Therefore for $`\delta ^2(t)=R^2(t)`$ one has $`\delta (t)\left(\frac{L_B}{\sqrt{3}}\frac{\sqrt{3}A}{2L_B}e^{t/\tau 1}\right)`$ The saturation value of $`\delta `$ is $`\delta _{max}=L_B/\sqrt{3}`$, so for $`t\tau _1`$, or equivalently for $`(\delta _{max}\delta )/\delta 1`$, we expect

$$\frac{d}{dt}\mathrm{ln}\delta =\lambda (\delta )=\frac{1}{\tau _1}\frac{\delta _{max}\delta }{\delta }$$
(B6)   
which is (12).

Let us remark that in the previous argument for $`\lambda (\delta )`$ for $`\delta \delta _{max}`$ the crucial point is the exponential relaxation to the asymptotic uniform distribution. In a generic deterministic chaotic system it is not possible to prove this property in a rigorous way. Nevertheless one can expect that this request is fulfilled at least in non-pathological cases. In chaotic systems the exponential relaxation to asymptotic distribution corresponds to have the second eigenvalue $`\alpha `$ of the Perron-Frobenius operator inside the unitary circle; the relaxation time is $`\tau _1=\mathrm{ln}|\alpha |`$ .

## C Synthetic turbulent velocity fields

The generation of a synthetic turbulent field which reproduces the relevant statistical features of fully developed turbulence is not an easy task. Indeed to obtain a physically sensible evolution for the velocity field one has to take into account the fact that each eddy is subject to the action of all other eddies. Actually the overall effect amounts only to two main contributions, namely the sweeping exerted by larger eddies and the shearing due to eddies of comparable size. This is indeed a substantial simplification, but nevertheless the problem of properly mimicking the effect of sweeping is still unsolved. 

To get rid of these difficulties we shall limit ourselves to the generation of a synthetic velocity field in Quasi-Lagrangian (QL) coordinates , thus moving to a frame of reference attached to a particle of fluid $`𝒓_1(t)`$. This choice bypasses the problem of sweeping, since it allows to work only with relative velocities, unaffected by advection. As a matter of fact there is a price to pay for the considerable advantage gained by discarding advection, and it is that only the problem of two-particle dispersion can be well managed within this framework. The properties of single-particle Lagrangian statistics cannot, on the contrary, be consistently treated.

The QL velocity differences are defined as

$$𝒗(𝒓,t)=𝒖(𝒓_1(t)+𝒓,t)𝒖(𝒓_1(t),t),$$
(C1)   
where the reference particle moves according to

$$\frac{\mathrm{d}𝒓_1(t)}{\mathrm{d}t}=𝒖(𝒓_1(t),t).$$
(C2)   
These velocity differences have the useful property that their single-time statistics are the same as the Eulerian ones whenever considering statistically stationary flows . For fully developed turbulent flows, in the inertial interval of length scales where both viscosity and forcing are negligible, the QL longitudinal velocity differences show the scaling behavior

$$\left|𝒗(𝒓)\frac{𝒓}{r}\right|^pr^{\zeta _p}$$
(C3)   
where the exponent $`\zeta _p`$ is a convex function of $`p`$, and $`\zeta _3=1`$. This scaling behavior is a distinctive statistical property of fully developed turbulent flows that we shall reproduce by means of a synthetic velocity field. In the QL reference frame the first particle is at rest in the origin and the second particle is at $`𝒓_2=𝒓_1+𝑹`$, advected with respect to the reference particle by the relative velocity

$$𝒗(𝑹,t)=𝒖(𝒓_1(t)+𝑹,t)𝒖(𝒓_1(t),t)$$
(C4)   
By this change of coordinates the problem of pair dispersion in an Eulerian velocity field has been reduced to the problem of single particle dispersion in the velocity difference field $`𝒗(𝒓,t)`$. This yields a substantial simplification: it is indeed sufficient to build a velocity difference field with proper scaling features in the radial direction only, that is along the line that joins the reference particle $`𝒓_1(t)`$ – at rest in the origin of the QL coordinates – to the second particle $`𝒓_2(t)=𝒓_1(t)+𝑹(t)`$. To appreciate this simplification, it must be noted that actually all moments of velocity differences $`𝒖(𝒓_1(t)+𝒓^{},t)𝒖(𝒓_1(t)+𝒓,t)=𝒗(𝒓^{},t)𝒗(𝒓,t)`$ should display power law scaling in $`|𝒓^{}𝒓|`$. Actually these latter differences never appear in the dynamics of pair separation, and so we can limit ourselves to fulfill the weaker request (C3). Needless to say, already for three particle dispersion one needs a field with proper scaling in all directions.

We limit ourselves to the two-dimensional case, where we can introduce a stream function for the QL velocity differences

$$𝒗(𝒓,t)=\times \psi (𝒓,t).$$
(C5)   
The extension to a three dimensional velocity field is not difficult but more expensive in terms of numerical resources.

Under isotropic conditions, the stream function can be decomposed in radial octaves as

$$\psi (𝒓,\theta ,t)=\underset{i=1}{\overset{N}{}}\underset{j=1}{\overset{n}{}}\frac{\varphi _{i,j}(t)}{k_i}F(k_ir)G_{i,j}(\theta ).$$
(C6)   
where $`k_i=2^i`$. Following a heuristic argument, one expects that at a given $`r`$ the stream function is essentially dominated by the contribution from the $`i`$ term such that $`r2^i`$. This locality of contributions suggests a simple choice for the functional dependencies of the “basis functions”:

$$F(x)=x^2(1x)\text{for}\mathrm{\hspace{0.17em}\hspace{0.17em}0}x1$$
(C7)   
and zero otherwise,

$$G_{i,1}(\theta )=1,G_{i,2}(\theta )=\mathrm{cos}(2\theta +\phi _i)$$
(C8)   
and $`G_{i,j}=0`$ for $`j>2`$ ($`\phi _i`$ is a quenched random phase). It is worth remarking that this choice is rather general because it can be derived from the lowest order expansion for small $`r`$ of a generic streamfunction in Quasi-Lagrangian coordinates.

It is easy to show that, under the usual locality conditions for infra red convergence, $`\zeta _p<p`$ , the leading contribution to the $`p`$-th order longitudinal structure function $`|v_r(r)|^p`$ stems from $`M`$-th term in the sum (C6), $`|v_r(r)|^p|\varphi _{M,2}|^p`$ with $`r2^M`$. If the $`\varphi _{i,j}(t)`$ are stochastic processes with characteristic times $`\tau _i=2^{2i/3}\tau _0`$, zero mean and $`|\varphi _{i,j}|^pk_i^{\zeta _p}`$, the scaling (C3) will be accomplished. An efficient way of to generate $`\varphi _{i,j}`$ is :

$$\varphi _{i,j}(t)=g_{i,j}(t)z_{1,j}(t)z_{2,j}(t)\mathrm{}z_{i,j}(t)$$
(C9)   
where the $`z_{k,j}`$ are independent, positive definite, identically distributed random processes with characteristic time $`\tau _k`$, while the $`g_{i,j}`$ are independent stochastic processes with zero mean, $`g_{i,j}^2k_i^{2/3}`$ and characteristic time $`\tau _i`$. The scaling exponents $`\zeta _p`$ are determined by the probability distribution of $`z_{i,j}`$ via

$$\zeta _p=\frac{p}{3}\mathrm{log}_2z^p.$$
(C10)   
As a last remark we note that by simply fixing the $`z_{i,j}=1`$ we recover the Kolmogorov scaling, which has been used in the simulations presented in section IV B

FIGURE CAPTIONS

a) Three realizations of $`R^2(t)`$ as a function of $`t`$ built as follows: $`R^2(t)=\delta _0^2\mathrm{exp}(2\gamma t)`$ if $`R^2(t)<1`$ and $`R^2(t)=2D(tt_{})`$ with $`\gamma =0.08,0.05,0.3`$ and $`\delta _0=10^7,D=1.5`$. b) $`R^2(t)`$ as function of $`t`$ averaged on the three realizations shown in figure 1a. The apparent anomalous regime and the diffusive one are shown. c) $`\lambda (\delta )`$ vs $`\delta `$, with Lyapunov and diffusive regimes.

Lagrangian motion given by the Rayleigh-Bénard convection model with: $`A=0.2,B=0.4,\omega =0.4,k=1.0,a=\pi `$, the number of realizations is $`𝒩=2000`$ and the series of thresholds is $`\delta _n=\delta _0r^n`$ with $`\delta _0=10^4`$ and $`r=1.05`$. $`\lambda (\delta )`$ vs $`\delta `$, in a closed domain with $`6`$ (crosses) and $`12`$ (diamonds) convective cells. The lines are respectively: (a) Lyapunov regime with $`\lambda =0.017`$; (b) diffusive regime with $`D=0.021`$; (c) saturation regime with $`\delta _{max}=19.7`$; (d) saturation regime with $`\delta _{max}=5.7`$.

(a)$`R^2(t)`$ for the four vortex system with $`\mathrm{\Gamma }_1=\mathrm{\Gamma }_2=\mathrm{\Gamma }_3=\mathrm{\Gamma }_4=1`$. The threshold parameter is $`r=1.03`$ and $`\delta _0=10^4`$, the dashed line is the power law $`R^2(t)t^{1.8}`$. The number of realizations is $`𝒩=2000`$. (b)$`\lambda (\delta )`$ vs $`\delta `$ for the same model and parameters. The horizontal line indicates the Lyapunov exponent ($`\lambda =0.14`$), the dashed curve is the saturation regime with $`\delta _{max}=0.76`$

FSLE computed for particle diffusion in a Sierpinsky Carpet of fractal dimension $`D_f=\mathrm{log}(8)/\mathrm{log}(3)`$ obtained by iteration of the unit structure up to a resolution $`3^5`$, one has: $`\lambda (\delta )\delta ^{1/.45}`$, which is in agreement with the value obtained for $`<R(t)>`$ versus $`t`$ shown in the inset (i.e. $`<R(t)>t^{.45}`$).

Relative dispersion of Lagrangian trajectories, $`<R_{x,y}^2(t)>`$ versus $`t`$, in the Adriatic Sea, for data (continuous line) and model (dashed line), along the natural axes of the basin: (a) transverse direction ($`x`$axis) and (b) longitudinal direction ($`y`$axis). The time is measured in $`days`$ and the mean square radius of the cluster is in $`km^2`$.

FSLE of Lagrangian trajectories in the Adriatic Sea, for data (continuous line) and model (dashed line). The scale $`\delta `$ is in $`km`$, $`\lambda (\delta )`$ is in $`days^1`$.

Relative dispersion $`R(t)`$ for $`N=20`$ octaves synthetic turbulent simulation averaged over $`10^4`$ realizations. The line is the theoretical Richardson scaling $`t^{3/2}`$.

Average inverse doubling time $`1/T(R)`$ for the same simulation of the previous figure. Observe the enhanced scaling region. The line is the theoretical Richardson scaling $`R^{2/3}`$.