Chaotic Lagrangian models for turbulent relative dispersion

TL;DR

This paper introduces a deterministic chaotic model for turbulent relative dispersion that avoids common stochastic simulation drawbacks, accurately analyzes dispersion statistics, and demonstrates efficiency across various Reynolds numbers.

Contribution

The paper presents a novel chaotic Lagrangian model for turbulence that removes the sweeping effect and accurately captures dispersion statistics using FSLE analysis.

Findings

Chaotic flows effectively model turbulent trajectories.

Model removes the sweeping effect via quasi-Lagrangian coordinates.

Numerical experiments show model efficiency across Reynolds numbers.

Abstract

A deterministic multi-scale dynamical system is introduced and discussed as prototype model for relative dispersion in stationary, homogeneous and isotropic turbulence. Unlike stochastic diffusion models, here trajectory transport and mixing properties are entirely controlled by Lagrangian Chaos. The anomalous "sweeping effect", a known drawback common to kinematic simulations, is removed thanks to the use of quasi-Lagrangian coordinates. Lagrangian dispersion statistics of the model are accurately analysed by computing the Finite-Scale Lyapunov Exponent (FSLE), which is the optimal measure of the scaling properties of dispersion. FSLE scaling exponents provide a severe test to decide whether model simulations are in agreement with theoretical expectations and/or observation. The results of our numerical experiments cover a wide range of "Reynolds numbers" and show that chaotic deterministic flows can be very efficient, and numerically low-cost, models of turbulent trajectories in stationary, homogeneous and isotropic conditions. The mathematics of the model is relatively simple and, in a geophysical context, potential applications may regard small-scale parametrization issues in general circulation models, mixed layer and/or boundary layer turbulence models as well as Lagrangian predictability studies.

Tables

Table 1. Table 1: Standard set-up for the 4-decade inertial range configuration analyzed in this paper. Parameters are expressed in non dimensional units.

	Kolmogorov constant	$C_K=1$
	turbulent dissipation rate	$ϵ={10}^{-3}$
	number of modes	$N_m=29$
	mode density	$a=\sqrt{2}$
	integral length	$L_0=2^{14}$
	wavelengths	$l_m=2L_0/a^{m-1}$ ($l_1=2L_0$)
	wave-numbers	$k_m=2\pi /l_m$
	square velocities	$A_m^{}{}_{}{}^2=2C_K{ϵ}^{2/3}k_m^{}{}_{}{}^{-5/3}\Delta k_m$ ($\Delta k_m=k_{m+1}-k_m$)
	turnover times	$t_m=2l_m/A_m$
	oscillation amplitudes	${\epsilon }_m=0.2l_m$
	oscillation pulsations	${\omega }_m=4\pi /t_m$
	Kolmogorov length	$\eta =l_{N_m}/2=1$
	integral time	$T_0=L_0^{2/3}{\left(2C_K\right)}^{-1/2}{ϵ}^{-1/3}$
	integration time step	$\Delta t={10}^{-3}{t}_{N_m}$

Equations

$$\lambda (\delta )\equiv \frac{1}{⟨\tau (\delta )⟩}\ln \frac{\varrho \cdot \delta }{\delta }.$$

$$\mathbf{u}=(\frac{\partial \Psi }{\partial x_2},-\frac{\partial \Psi }{\partial x_1})$$

$$\frac{d\mathbf{r}}{dt}=\mathbf{u}$$

$$⟨||\delta \mathbf{x}(t)||⟩\sim ||\delta \mathbf{x}(0)||{\mathrm{e}}^{{\lambda }_{L}t}$$

$$u_{i_1}=\frac{\partial {\Psi }_{i_3}}{\partial x_{i_2}}-\frac{\partial {\Psi }_{i_2}}{\partial x_{i_3}}$$

$$s_i(t)\equiv x_i-\epsilon \sin ({\omega }_{i}t)$$

$${\Psi }_{i_1}(x_{i_2},x_{i_3},t)=\frac{A}{k}\sin [k{s}_{i_2}(t)]\sin [k{s}_{i_3}(t)]$$

$$u_{i_1}[\mathbf{x},t]=A\{\sin [k{s}_{i_1}(t)]\cos [k{s}_{i_2}(t)]-\sin [k{s}_{i_1}(t)]\cos [k{s}_{i_3}(t)]\}$$

$$⟨||\delta \mathbf{x}(t)||⟩\sim ||\delta \mathbf{x}(0)||{\mathrm{e}}^{{\lambda }_{L}t},$$

$$⟨||\delta \mathbf{x}(t)|{|}^2⟩\sim 4D_{l}t,$$

$${\Psi }_{i_1}(x_{i_2},x_{i_3},t)=\sum _{m=1}^{N_m}{\Psi }_{i_1}^{(m)}(x_{i_2},x_{i_3},t)$$

$$A_m^2\equiv 2C_K{ϵ}^{2/3}{k}_m^{-5/3}\Delta k_m$$

$${\mathbf{x}}^{(C)}\equiv \frac{{\mathbf{x}}^{(1)}+{\mathbf{x}}^{(2)}}{2}$$

$$u_{i_1}^{(m)}[{\mathbf{x}}^{(j)},t]=A_m\{\sin [k_m{s}_{i_1}^{(j)}(t)]\cos [k_m{s}_{i_2}^{(j)}(t)]-\sin [k_m{s}_{i_1}^{(j)}(t)]\cos [k_m{s}_{i_3}^{(j)}(t)]\}$$

$${\mathbf{s}}^{(j)}(t)\equiv {\mathbf{x}}^{(j)}-\mathbf{\epsilon }\sin (\omega t)$$

$${\mathbf{s}}_{QL}^{(1,2)}(t)\equiv [{\mathbf{x}}^{(1,2)}-{\mathbf{x}}^{(C)}]-\mathbf{\epsilon }\sin (\omega t)$$

$$\{\begin{array}{l}\mathbf{u}[{\mathbf{x}}^{(1,2)},t]={\mathbf{u}}^{(1)}[{\mathbf{s}}^{(1,2)}(t)]+\sum _{m=2}^{N_m}{\mathbf{u}}^{(m)}[{\mathbf{s}}_{QL}^{(1,2)}(t)]\\ \\ \frac{d{\mathbf{x}}^{(1,2)}}{dt}=\mathbf{u}[{\mathbf{x}}^{(1,2)},t]\end{array}$$

$$p_R^{\ast }(r,t)\simeq \frac{C}{C_Rϵt^3}\exp \left(-C^{\prime }\frac{r^{2/3}}{C_R^{1/3}{ϵ}^{1/3}t}\right)$$

$$p_e^{\ast }(\tau )\simeq \exp \left(-C^{\prime \prime }\frac{\tau }{⟨\tau (\delta )⟩}\right)$$

Full text

Kinematic Models, Deterministic chaos, Lagrangian Turbulence, Finite-Scale Lyapunov Exponent

I Introduction

The physical characteristics of a fluid dynamical system can be classified in two major categories: Eulerian, dealing with vector and scalar fields as functions of space and time coordinates, and Lagrangian, dealing with quantities related to the motion of fluid particles. As far as transport processes are concerned, the Lagrangian approach provides information not trivially obtainable from the only knowledge of the velocity field (Crisanti et al., 1991). The existence of Lagrangian Chaos (Aref, 1984; Ottino, 1989), has shown the possibility to have efficient transport and dispersion even in regular velocity fields, e.g. periodic in space and/or time. In addition to this, turbulence (Frisch, 1995) can act further as mechanism of particle dispersion on intermediate scales of motion.

Modelling Lagrangian particle trajectory evolution is, therefore, an intriguing task both from the theoretical and methodological points of view. Ideally, starting from the knowledge of the velocity field, in every point and at every time, one can numerically integrate fluid particle trajectories by means of a computer over arbitrarily long time intervals. Of course, such an approach cannot be used in realistic situations, since models and/or observative data are always affected by finite resolution. Trying to accurately describe the dynamics at all scales of motion presents, in general, insurmountable difficulties, even in the hypothesis of having a full knowledge of the physics of the system (which, in actual fact, is only partly true). So, usually, one has to find a suitable compromise between the level of resolution of the model and the number of degrees of freedom to take into account.

As far as Lagrangian dynamics is concerned, i.e. the evolution of fluid particle trajectories, there are two major statistics to consider: one-particle, or absolute, dispersion, i.e. the mean displacement from the initial positions, which depends essentially on the large-scale dynamical features, and two-particle, or relative, dispersion, i.e. the mean separation between trajectories, which carries on interesting information about the physics of the system (LaCasce, 2008). In the present paper we will deal only with relative dispersion modelling. At this regard, we recall that, in some sense one could say luckily, the Lagrangian properties of a fluid, e.g. mixing and diffusion, do not depend strongly on the details of the velocity field, but do depend, indeed, mainly on the relationship between the characteristic spatial and temporal scales of the system (Crisanti et al., 1991). As an example, if we consider the diffusion process simulated by a $\tau-$correlated Langevin equation (as in the classical Brownian motion model), such that $\sigma^{2}\tau\sim D$, where $D$ is the diffusion coefficient, $\sigma^{2}$ is the variance of the velocity fluctuations and $\tau$ is the velocity auto-correlation time, we notice that the same type of dispersion process, i.e. ballistic regime for times smaller than $\tau$ and asymptotic diffusive regime for times larger than $\tau$, can be simulated by the chaotic, deterministic scattering across a lattice of (unsteady) kinematic eddies, as long as $l^{2}/\tau\sim D$, where $l$ and $\tau$ indicate the eddy size and turnover time, respectively. We will come back to this point later. Such a fact opens interesting possibilities from the Lagrangian modelling viewpoint, since there is no need to compute a realistic turbulent flow if the scope is only that to simulate the the relative dispersion process in turbulent flows. For example, synthetic or kinematic models can accomplish this task with good efficiency and low computational cost, provided some important aspects are suitable considered. The scope of this work is to present and discuss a general kinematic model that can be utilized to numerically simulate Lagrangian trajectories in an ideal turbulent environment in which stationarity, isotropy and homogeneity conditions are assumed to be fulfilled. What is known as dynamical system approach to Lagrangian transport and mixing (Aref, 1984; Ottino, 1989; Crisanti et al., 1991; Bohr et al., 1998; Ott, 2002; Weiss and Provenzale, 2008) will be adopted as guiding line throughout this paper.

The idea of implementing kinematic simulations of turbulence, as alternative to stochastic Lagrangian modelling based on Langevin equations (Thomson, 1987), is not new, and many papers can be counted in support of this strategy (Fung et al., 1992; Fung and Vassilicos, 1998; Nicolleau and Vassilicos, 2003; Osborne et al., 2006; Nicolleau and Nowakowski, 2011). Although our kinematic model differs from similar models used in the cited literature, since it is purely deterministic, the philosophy of the approach is quite the same: using an analytical velocity field to numerically compute trajectories that behave, in some statistical sense, like a real turbulent flow. Some Authors (Thomson and Devenish, 2005; Devenish and Thomson, 2009; Eyink and Benveniste, 2013) have raised objections as far as the use of kinematic simulations is concerned. It was correctly argued, indeed, that a kinematic velocity field, characterized by structures (eddies) that do not drift all over the domain, therefore opposing an unnatural resistance to the dragging force exerted by the large-scale energetic field, cannot generate trajectories that simulate turbulent diffusion in the correct way. It is known that pair dispersion, within the inertial range of scales where the turbulent cascade has developed, is expected to evolve according to the so-called locality hypothesis, i.e. it is assumed that only those structures having characteristic size of the same order as the particle separation scale contribute efficiently to the dispersion process (Frisch, 1995). This means that, e.g., the mean square particle displacement is expected to obey the well known ”$t^{3}$” Richardson-Obukhov law (Richardson, 1926).

In a kinematic model, this property can be obscured (or distorted) by the so-called sweeping effect due to the most energetic components of the velocity field. These large-scale structures drag the particle pairs across the domain but do not have any effect on the kinematic eddies at smaller scales. This implies that a particle pair does not spend sufficient time around a local eddy such as to “thermalize” to the local dynamics. As a consequence, the mean square relative dispersion deviates from the “$t^{3}$” law of an amount that grows macroscopically in the limit of long inertial range and high time resolution, as explained and discussed in detail by Thomson and Devenish (2005) (Thomson and Devenish, 2005).

In the meanwhile, some years ago, Lacorata et al. (Lacorata et al., 2008) have first applied a type of ”correct” kinematic model as sub-grid-scale parameterization of two-particle dispersion in a Large-Eddy Simulation of planetary boundary layer turbulence. Later, this methodology has been extended also to oceanographic applications (Palatella et al., 2014; Lacorata et al., 2014). The aim of the present work is to offer a further and specific contribution to clarify this question. A very general formulation of a three-dimensional, multi-scale, deterministic kinematic model will be introduced and analyzed in detail by means of numerical simulations. The Lagrangian turbulence model is built up assuming the validity of K41 theory (Frisch, 1995). For many applications, this does not have to be considered as a limitation. A severe test on the skills of the model is performed by computing the Finite-Scale Lyapunov Exponents (see next section), which measure the spectrum of relative dispersion rates at all scales of motion and allow a direct comparison with the theory.

This paper is organized as follows: the FSLE technique is recalled in Section II; in Section III, the Kinematic Lagrangian Model is introduced and discussed; the results of the numerical simulations are described in Section IV and a summary of the conclusions that can be drawn at the end of this work is reported in Section V.

II Relative dispersion: the Finite-Scale Lyapunov Exponent

The Finite-Scale Lyapunov Exponent (FSLE) is a scale-dependent measure of the mean growth rate of the distance between two trajectories, an established technique employed in a wide range of applications (Aurell et al., 1996, 1997; Boffetta and Celani, 2000; LaCasce and Ohlmann, 2003; Lacorata et al., 2004; LaCasce, 2008; Berti et al., 2011; Lacorata and Espa, 2012; Espa et al., 2014), from studies on the growth of finite perturbations in chaotic dynamical systems to turbulent pair dispersion in ocean and atmosphere, as well as in fluid dynamics laboratory experiments. The idea of measuring the diffusivity, i.e. the time derivative of the mean square particle displacement, as function of the particle separation actually dates back to the times of Richardson (Richardson, 1926). The FSLE adopts the same philosophy, i.e. it measures the dispersion rate at fixed scale, thus avoiding issues related to the fixed time averaging procedure (Boffetta et al., 2000). Although definition and properties of the FSLE have already been described and discussed in detail in literature, e.g., see (Cencini and Vulpiani, 2013) for a very general review, for the sake of self-consistency we recall here its basic characteristics. Given two trajectories, separated at time $t$ by a distance $r(t)$, let us define $\tau$ as the time interval during which the distance grows from $\delta$ to $\varrho\cdot\delta$, where $\varrho\gtrsim 1$. Keeping $\delta$ and $\varrho$ fixed, the mean growth time $\langle\tau\rangle$ from $\delta$ to $\varrho\cdot\delta$, suitably averaged (Cencini and Vulpiani, 2013), defines the FSLE $\lambda(\delta)$ according to the formula:

[TABLE]

The quantity $\tau(\delta)$ is the first exit time of the distance $r$ from scale $\delta$ to scale $\varrho\cdot\delta$. It is implicitly assumed that $\langle\tau(\delta)\rangle$, i.e. the phase space average or, equivalently, the average over an arbitrarily large number of numerical experiments, is a function of $\delta$, for $\varrho\sim O(1)$, but it does not depend on $t$, under the hypothesis of stationary statistics. The idea is to define $N_{s}+1$ scales, such that $\delta_{n}=\varrho\cdot\delta_{n-1}$, with $n=1,...,N_{s}$, and to compute the FSLE (1) for each of them. In a Lagrangian dynamical context, the phase space coincides with the familiar physical space and the trajectories are those of passive tracer particles in a given velocity field. Normally, the initial scale $\delta_{0}$ should be smaller than the least characteristic length of the flow, e.g. the Kolmogorov scale, and the final scale $\delta_{N_{s}}$ should be of the order of the integral scale of the flow or possibly larger, depending on the size of the system. A common value of the amplification ratio, adopted also in this work, is $\varrho=\sqrt{2}$. It is worth stressing a couple of remarks: the time $\tau(\delta)$ is to be computed as the first crossing time of the scale $\varrho\cdot\delta$, starting from the scale $\delta$; the factor $\varrho$ must be chosen not too close to unity, in order to avoid saturation issues of the FSLE at small scale separations due to the finite time resolution of the trajectories111The algorithm used for the computation of the FSLE is anyway prepared to compensate possible clipping effects (Boffetta et al., 2000), and not too larger than unity, otherwise an accurate scale dependent description of the dispersion regimes could not be feasible.

For a fully developed turbulent flow, in which a self-similar energy cascade establishes in the inertial sub-range between a forcing scale $L_{0}$ (integral length scale) and a dissipative scale $\eta$ (Kolmogorov length scale), three major regimes are expected to be observable:

a) exponential separation, $\langle r(t)^{2}\rangle\sim r(0)^{2}{\rm e}^{2\lambda_{L}t}$, on scales smaller than the Kolmogorov length, $r\ll\eta$, where $\lambda_{L}$ is the Lagrangian Maximum Lyapunov Exponent222Taking into account intermittent fluctuations of the effective Lyapunov exponent one has $\langle r(t)^{2}\rangle\sim r(0)^{2}{\rm e}^{L(2)\,t}$, where $L(2)\geq 2\lambda_{L}$. In the present work we will ignore this detail. of the flow (Crisanti et al., 1991);

b) Richardson-Obukhov scaling $\langle r(t)^{2}\rangle\sim C_{R}\,\epsilon\,t^{3}$ in the inertial range, $\eta<r<L_{0}$, where $\epsilon$ is the rate of energy dissipation of the cascade and $C_{R}$ is the non-dimensional Richardson’s constant (Richardson, 1926); henceforth we will refer to any power law of the type $\langle r(t)^{2}\rangle\sim t^{\nu}$, for which $\nu>1$, as super-diffusion or, equivalently, super-diffusive regime;

c) Standard diffusion $\langle r(t)^{2}\rangle\sim 4\,D_{E}\,t$ on scales larger than the integral length, $r\gg L_{0}$, where $D_{E}$ is the eddy diffusion coefficient (Taylor, 1921).

Let us briefly discuss each of these regimes and the relationship with the scaling properties of the FSLE, see Figure 1:

a) the small-scale exponential separation is characterized by a scale-independent mean growth rate of the distance between two trajectories. In terms of FSLE this corresponds to $\lambda(\delta)=\lambda_{L}$ for $\delta\ll\eta$;

b) by dimensional arguments it can be shown that, inside the inertial range $\eta<\delta<L_{0}$, the Richardson-Obukhov law corresponds to the scaling $\lambda(\delta)\sim(C_{R}\cdot\epsilon)^{1/3}\,\delta^{-2/3}$;

c) as for the previous point, by dimensional arguments it is possible to show that, in the limit of large separation scales $\delta\gg L_{0}$, standard diffusion corresponds to the scaling $\lambda(\delta)\sim D_{E}\,\delta^{-2}$.

From the FSLE analysis it is possible not only to distinguish the type of regime a), b) or c) in correspondence of a given range of scales, but also to estimate, at least as order of magnitude, the characteristic parameters of each regime, $\lambda_{L}$, $\epsilon$ or $D_{E}$. At this regard, it is clear that an unequivocal measure of the scaling exponents, characteristic of the relative dispersion process, is needed in order to assess whether or not numerical experiments are in agreement with the theory. We would like to stress again that the FSLE approach, where the ”proper” variable is $\delta$, is rather different than looking at $\langle r^{2}(t)\rangle$ vs $t$. In the latter case often, at a given time, one can have contamination of effects due to different scales. As a result of such a contamination, $\langle r^{2}(t)\rangle$ can strongly depend on $\langle r^{2}(0)\rangle$ (Cencini and Vulpiani, 2013). This drawback is ruled out by the FSLE analysis technique.

III Kinematic Lagrangian Model

III.1 Velocity field

Let us assume that the flow is incompressible. For a 2D velocity field, ${\bf u}=\{v_{1},v_{2}\}$ is a continuous function of the spatial coordinates ${\bf x}=\{x_{1},x_{2}\}$ and of the time $t$. The incompressibility condition $\nabla\cdot{\bf u}=0$ can be fulfilled by introducing the stream-function $\Psi(x_{1},x_{2},t)$ and, then, writing the velocity field as:

[TABLE]

The evolution of the Lagrangian coordinates of a fluid particle is given by:

[TABLE]

The dynamical system (3), with the velocity field defined in (2), is formally a Hamiltonian system with $\Psi$ playing the role of the Hamiltonian function (Ottino, 1989). Once the stream-function, or equivalently the velocity field, is assigned, the system formed by (2) and (3) is called kinematic model. Lagrangian trajectories computed in a 2D kinematic model are regular, i.e. non chaotic, if ${\bf u}$ does not depend explicitly on time $t$. If, conversely, ${\bf u}$ is a nonlinear vector field explicitly depending on the time $t$, Lagrangian dynamics (3) can be chaotic, i.e. two arbitrarily close trajectories tend to separate exponentially in time:

[TABLE]

in the limit $||\delta{\bf x}||\to 0$, where $||\delta{\bf x}||\equiv||{\bf x}^{(1)}-{\bf x}^{(2)}||$ is the distance between two trajectories ${\bf x}^{(1)}$, ${\bf x}^{(2)}$, and $\lambda_{L}>0$ is the Maximum Lagrangian Lyapunov Exponent (LLE) of the flow. The LLE $\lambda_{L}$ can be positive even when ${\bf u}$ is a regular function, e.g. periodic in space and time. A paradigmatic example of a 2D chaotic flow is the Rayleigh-Bénard convection model by Solomon and Gollub (Solomon and Gollub, 1988; Crisanti et al., 1991).

When perturbing an integrable Hamiltonian system, the transition from non chaotic to full chaotic regime occurs via progressive destruction of the KAM tori, until a fully chaotic regime is reached for suitable values of the perturbation parameters, i.e. in the so-called overlap of the resonances regime (Chirikov, 1979). Therefore, typically, chaotic motions co-exist with islands of regular motions in the phase space, and the maximum Lyapunov exponent is not positive everywhere but, averaging over all possible initial conditions, it has positive mean (Ottino, 1989; Ott, 2002). In the case discussed in the present work, the model velocity field is suitably perturbed in order to reach the Chirikov regime, without trapping regions or other kind of regular islands.

Let us consider, now, the general three dimensional case. Let ${\bf\Psi}=\{\Psi_{1},\Psi_{2},\Psi_{3}\}$ be a continuous and differentiable vector field, function of the spatial coordinates, ${\bf x}=\{x_{1},x_{2},x_{3}\}$, and of the time, $t$. Let ${\bf u}=\{u_{1},u_{2},u_{3}\}$ be a velocity field defined as ${\bf u}={\bf\nabla}\times{\bf\Psi}$, i.e.:

[TABLE]

where $\{i_{1},i_{2},i_{3}\}$ stand for all cyclic permutations of $\{1,2,3\}$. Of course, ${\bf u}$ is a solenoidal field, i.e. ${\bf\nabla}\cdot{\bf u}\equiv 0$. It must be observed that, in the three dimensional case, providing ${\bf u}$ is non linear, a dynamical system of the type (3) can display Lagrangian chaos even when ${\bf u}$ is not explicitly time dependent (Crisanti et al., 1991). Henceforth we will consider only non linear velocity fields that are periodic functions both in space and time. At this regard, for simplicity of notation, let us define the following change of variables from ${\bf x}=\{x_{1},x_{2},x_{3}\}$ to ${\bf s}(t)=\{s_{1}(t),s_{2}(t),s_{3}(t)\}$ as:

[TABLE]

where $i=\{1,2,3\}$, ${\bf x}$ is the position vector in the fixed reference frame, and ${\bf s}(t)$ oscillates periodically around ${\bf x}$ with frequency ${\bf\omega}/(2\pi)$ and amplitude $\varepsilon$. For technical reasons that will be clear later it is convenient to define the same amplitude $\varepsilon$ for all components, and the pulsations $\omega_{1}\simeq\omega_{2}\simeq\omega_{3}\simeq\omega$ nearly equal to one another. Let us now define a 3D periodic lattice of convective cells as:

[TABLE]

where, as for Eq. (5), the set $\{i_{1},i_{2},i_{3}\}$ is a cyclic permutation of $\{1,2,3\}$; $k=2\pi/l$ and $l$ are spatial wave-number and wavelength, respectively; $\varepsilon$ and $\omega/(2\pi)$, contained in the time dependent terms, $s_{1}(t)$, $s_{2}(t)$ and $s_{3}(t)$, are amplitude and frequency of the stream-function oscillation. The parameter $A$ defines the velocity scale of the flow. From (5) and (7) the three components of the kinematic velocity field, adopting the cyclic index notation as above, turn out to be:

[TABLE]

The field (8) can be seen as a three dimensional version of the Solomon and Gollub model (Solomon and Gollub, 1988). Other options are possible, of course, in the choice of the analytical form of ${\bf u}$. The ABC flow (Crisanti et al., 1991), for example, or the so-called Double Stream Function (DSF) model (Lacorata et al., 2008) are other good candidates. In analogy with the overlap of the resonances in Hamiltonian dynamics (Chirikov, 1979), we must expect that, for a certain range of values of the lattice oscillation parameters $\varepsilon$ and $\omega$, every possible bounded region of regular motion is destroyed and Lagrangian trajectories can evolve chaotically across the whole domain. In other words, for any $\omega$, $\varepsilon$ must be larger than a certain critical threshold $\varepsilon_{c}(\omega)$ in order to let the model reach a suitable working point, i.e. a regime in which Lagrangian chaos attains a good efficiency as mechanism of trajectory dispersion. At this regard, the standard set up of the parameters can be defined on the basis of the following considerations. The 3D periodic lattice is formed by cubic elementary cells of edge $l/2$, and characteristic time scale $t_{c}=2\,l\,/\,A$. The periodic oscillations of the structure let a fluid particle escape from its initial cell and perform a sort of random walk through the lattice, with characteristic time of the same order as $t_{c}$, providing the oscillation frequency is of the same order as the turnover frequency, i.e. $\omega/(2\pi)\sim 1/t_{c}$. The oscillation amplitude can be fixed to some fraction of the characteristic length of the cells, $\varepsilon/l\sim 10^{-1}$.

In all cellular flows of this type, characterized by only one scale of motion, relative dispersion displays two regimes: small-scale exponential separation,

[TABLE]

in the limit $||\delta{\bf x}||\ll l$, with $\lambda_{L}\sim O(1/t_{c})$, and large-scale standard diffusion,

[TABLE]

for $||\delta{\bf x}||\gg l$ (Crisanti et al., 1991). In order to simulate an intermediate regime of turbulent dispersion, between these two asymptotic limits, it is mandatory to include different scales of motions in the velocity field. This can be attained by the superposition of a series of $N_{m}$ self-similar spatial modes, each corresponding to a different wave-number. With the cyclic index notation the resulting multi-scale stream-function is:

[TABLE]

In (9) the terms in the right-hand side have the same form as in (7) but now the parameters depend on the $m$ mode: $A_{m}$, $k_{m}$, $l_{m}$, $\varepsilon_{m}$ and $\omega_{m}$. The three components of the kinematic velocity field are now obtained from (5) and (9). At this point, the 3D kinematic model contains $N_{m}$ spatial modes, each having wavelength $l_{m}$ and characteristic velocity $A_{m}$. The wavelengths can be related to each other by a recursive rule: $l_{m}=L_{0}/a^{(m-1)}$, for $m=1,...,N_{m}$, and $a>1$. We let the largest scale $L_{0}$ correspond to the integral length and the smallest one $l_{N_{m}}\equiv\eta$ to the Kolmogorov length of the flow. With $L_{0}$ and $\eta$ fixed, the parameter $a$ determines the density of the modes. In order to simulate 3D isotropic homogeneous turbulence, the further step is to assign the Kolmogorov scaling (Bohr et al., 1998) to the velocities of the modes:

[TABLE]

where $k_{m}\equiv 2\pi/l_{m}$, $\Delta k_{m}\equiv k_{m+1}-k_{m}$, $m=1,...,N_{m}$; $\epsilon$ is the mean rate of turbulent dissipation and $C_{K}$ is a non dimensional adjustable constant (the equivalent Kolmogorov constant). Eq. (10) determines how energy is distributed among the spatial modes in the inertial range. Amplitude and frequency of the time dependent oscillating terms are set according to $\varepsilon_{m}/l_{m}\sim O(10^{-1})$ and $\omega_{m}t_{m}/(2\pi)\sim O(1)$, where the “local” time scales are defined as $t_{m}\equiv 2l_{m}/A_{m}$. The values of $\omega_{m}$ along the three directions $x_{1}$, $x_{2}$ and $x_{3}$, differ a little from each other of “irrational” factors numerically very close to 1. This is to avoid possible (even though unlikely) anomalous trapping of particles inside a box which might affect the chaotic diffusive process. At this point we have defined a 3D unsteady lattice of convective cells of various wavelengths ranging from $\eta$ to $L_{0}$, in which the velocities of spatial modes are related to the corresponding characteristic lengths by the Kolmogorov scaling. Integrating Lagrangian trajectories with this kinematic model and looking at the two-particle dispersion, we observe that, despite the scaling (10), the mean square separation, or equivalently the FSLE, departs from the expected Richardson-Obukhov scaling inside the inertial range, with a discrepancy increasing with the width of the inertial range, $L_{0}/\eta\gg 1$. As argued by Thomson and Devenish (Thomson and Devenish, 2005), this is due to the fact that the large-scale velocity field advects a particle pair, having an initially small separation, through the domain without advecting the surrounding small-scale structures, the so-called “sweeping” problem. As direct consequence, not only the locality hypothesis inside the inertial range is violated, but also the mean rate of exponential separation, at scales smaller than the Kolmogorov length, is sensitive to the integral scale and, ultimately, depends on the width of the inertial range.

The Maximum Lyapunov Exponent, on the contrary, is expected to be related to the inverse characteristic time of the smallest coherent structure of the flow (Bohr et al., 1998), i.e. $\lambda_{L}\sim 1/t_{\eta}$ where $t_{\eta}\sim\epsilon^{-1/3}\eta^{2/3}$ is the turnover time at the Kolmogorov scale $\eta$. The locality assumption, also, as far as the inertial range is concerned, is disregarded since particle pairs cannot “thermalize” to the local dynamics spending a sufficient time around eddies of the same size as the separation distance. The result is that the LLE depends on the large-scale velocity, with all other parameters fixed, and relative dispersion within the inertial range follows an anomalous scaling $\langle\delta{\bf x}(t)^{2}\rangle\sim t^{\nu}$ with $9/2\leq\nu\leq 6$, or, equivalently, in terms of FSLE $\lambda(\delta)\sim\delta^{-\mu}$ with $1/3\leq\mu\leq 4/9$. The reasons why the scaling exponents are expected to vary between these limits have been explained by Thomson and Devenish (Thomson and Devenish, 2005) and will not be repeated here. We are only concerned in evaluating under what conditions inertial range relative dispersion of the model is in agreement with Richardson’s law or not.

III.2 Quasi-Lagrangian coordinates

In order to treat the “sweeping” problem in a proper way, all spatial modes except the largest one can be written as functions of two particle (relative), instead of single particle (absolute), coordinates, i.e. making use of the so-called Quasi-Lagrangian frame technique (Boffetta et al., 1999; Lacorata et al., 2008). The point is to consider two particles at the same time and advect them together. The trajectories of a whole set of tracer particles can be, in this way, integrated pair by pair.

This strategy allows to fix the problem of having flow structures that do not move together with the fluid particles: computing the kinematic velocity field in the reference frame anchored to the center of mass of two particles simulates the simultaneous advection of kinematic eddies and particle pair. This advection may generated either by an external large-scale velocity field or, simply, by letting the largest velocity mode of the model be a function of single particle coordinates, as is the case treated in this work.

Let ${\bf x}^{(1)}$ and ${\bf x}^{(2)}$ be two fluid particle positions in the fixed reference frame. The vector position of their center of mass is:

[TABLE]

The $m$ mode velocity component in ${\bf x}^{(j)}$ ($j=1,2$) can be written as:

[TABLE]

where the usual cyclic index notation is adopted, and

[TABLE]

with ${\boldsymbol{\varepsilon}}\equiv\varepsilon(1,1,1)$. Let us consider both particles of a pair, ${\bf x}^{(1,2)}$, and rewrite (13) in the reference frame of their center of mass:

[TABLE]

where ${\bf x}^{(C)}$ is defined in (11). Let us write the full kinematic Lagrangian model as follows:

[TABLE]

We will call (15) QL-model, in which every mode except one ($m=1$) is written as function of the Quasi-Lagrangian coordinates relatively to a given particle pair. In other terms, the QL-model velocity field depends on the coordinates of both particles of a pair, via their center of mass. The kinematic field in (15) for ${\bf s}_{QL}\equiv{\bf s}$, i.e. setting ${\bf x}_{C}\equiv 0$ in (14), will be called E-model, in which every mode is written in terms of single particle coordinates ${\bf x}$. It is worth stressing that the formulation in terms of two particle coordinates does not imply that the QL-model is only a two-particle model (Boffetta et al., 1999), since the most energetic mode ($m=1$) in (15), function of single particle coordinates, simulates a large-scale advection acting on the single trajectories. We would like to remark that, in many applications, the role of large-scale advection is normally played, for example, by a general circulation model velocity field, so that there is no need to keep the $m=1$ mode as function of single particle coordinates any longer. In the next section, the results obtained from the Lagrangian simulations performed with the QL-model and the E-model will show the differences between the two types of configuration, as far as relative dispersion is concerned.

IV Results

A standard parameter set up of the 4-decade kinematic model is reported in Tab. 1. The width of the inertial range, in all analyzed cases, is determined by the integral length $L_{0}$, at fixed Kolmogorov length $\eta$. All other parameters are kept fixed in all numerical simulations.

We first checked the properties of relative dispersion in the time domain for the QL-model. We computed the mean square pair separation $\langle r(t)^{2}\rangle$, for the 4-decade configuration, to assess: the existence of the Richardson’s scaling law inside the inertial range, at given $\epsilon$; and an empirical relationship between the Richardson’s constant $C_{R}$ and the Kolmogorov constant $C_{K}$ of the model. These results are reported in Fig. 2. We observe that: the QL-model reproduces the expected Richardson’s regime for relative dispersion inside the inertial range, at various $C_{K}$ values; Richardson’s constant $C_{R}$ can be measured and reported as function of $C_{K}$, at given $\epsilon$; $C_{R}$ attains values in the expected range, $[0.1-1]$, for $C_{K}$ close to the theoretical value, $C_{K}\simeq 1$. Therefore, in all simulations we fixed $C_{K}=1$. We remark also that the empirical scaling relation between the two constants, $C_{R}\sim C_{K}^{3/2}$, can be justified on the basis of theoretical arguments by noticing that both quantities, $C_{R}$ and $C_{K}^{3/2}$, appear as rescaling factor of $\epsilon$ in the Richardson’s law and in the Kolmogorov spectrum, respectively.

FSLE’s for E-model and QL-model are plotted in Fig. 3, for identical parameter set up, varying the width of the inertial range over four orders of magnitudes. All curves, normalized to $\lambda_{L}$, collapse perfectly in the QL-model, fulfilling the invariance property of the Richardson scaling with respect to the inertial range. On the contrary, in the E-model such a property does not hold and the scaling deviates progressively from the Richardson’s law at increasing $L_{0}/\eta$ ratio. Outside the inertial range, at scales larger than the integral length, the scaling approaches a standard diffusive regime characterized by the same eddy diffusion coefficient, of course $L_{0}-$dependent, for both models. At small scales, i.e. below the Kolmogorov length, exponential separation appears but with a mean rate that, in the QL-model, is constant in all numerical experiments while, in the E-model, is sensitive to the integral scale. This picture confirms unequivocally that the conjecture advanced by Thomson and Devenish (2005) is verified for what concerns the behavior of kinematic simulations without the ”sweeping effect” correction, and, at the same time, that the Quasi-Lagrangian coordinate technique works perfectly to restore the right Richardson’s law on arbitrarily long inertial range.

A further evidence is shown in Fig. 4 where, in order to best highlight the differences, the FSLE’s are computed for an eight decade long inertial range in both models.

The numerical simulations have been performed with both fixed and adaptive integration time step $\Delta t$, showing no difference between the two methods. In the fixed time step case, $\Delta t$ is set to a value much smaller than the Kolmogorov time $t(\eta)$. This choice is very conservative since it permits a time resolution of the dispersion process that grows progressively with the scale of motion but, of course, it is also very expensive in terms of computation time when the inertial range attains several decades. On the other hand, the adaptive time step, defined as $\Delta t_{n}=10^{-2}\,l_{n}/A_{1}$ for $l_{n+1}<r<l_{n}$, is more efficient in terms of resource consuming and give exactly the same results of the former case, regardless the width of the inertial range. We remark that $l_{n}/A_{1}$ is the advecting time at scale $l_{n}$, and the condition for $dt$ is to be sufficiently small as to resolve well the dynamics of dispersion at local scale.

Last, we have studied the statistics of the dispersion process within the inertial range of the model (with the $L_{0}/\eta\sim 10^{4}$ set up, see Fig. 3) and checked it out with the theoretical predictions. At this regard, let us define the probability distribution function (pdf) of the distance $r$ between two particles at time $t$ as $p_{R}(r,t)$, and the pdf of the exit time $\tau$ at scale $\delta$ as $p_{e}(\tau)$. Following Boffetta and Celani (2000), we recall that, in the Richardson regime, the expected forms for the two pdf’s are, respectively:

[TABLE]

[TABLE]

where $\epsilon$ is the mean turbulent dissipation rate, $C_{R}$ is the Richardson’s constant, $\langle\tau(\delta)\rangle$ is the scale-dependent mean exit time, $C$, $C^{\prime}$ and $C^{\prime\prime}$ are constants whose precise value is not particularly important in this context.

The results are shown in Fig. 5. The pdf $p_{R}(r,t)$ is computed at three different times after the release of the particle pairs and, except for the core of the distribution around $r\simeq 0$, the tail is well fitted by the Richardson’s prediction (16) up to 4 or 5 times the standard deviation. The pdf $p_{e}(\tau)$ is computed for three different separation thresholds and, again, except for a limited range of $\tau\simeq 0$, the exponential tail is in agreement with the theoretical prediction (17) as confirmed by the collapse of the re-normalized curves (Boffetta and Celani, 2000). The discrepancy observed at small $r$ and $\tau$ is likely due to the ”smoothness” of the kinematic model at all scale separations with respect to the ”roughness” of a real turbulent field.

V Discussion and conclusions

In this paper, the problem of modelling two-particle dispersion in a turbulent fluid has been addressed by means of a dynamical system approach. Unlike other types of kinematic simulations, Lagrangian chaos and deterministic dynamics are the key elements characterizing the model we have here introduced and analyzed. The main conclusions at the end of this work can be summarized as follows:

1. Multi-scale dynamics typical of turbulent flows is guaranteed by the superposition of a number of self-similar modes.

2. Kolmogorov relationship between space and velocity parameters, scale by scale, define an inertial range in the model.

3. Quasi-Lagrangian coordinates assure the statistics of turbulent relative dispersion turn out to be in very good agreement with the theoretical expectations.

While points 1) and 2) are already established facts in literature, point 3) is an original element, fundamental to the correct behavior of the kinematic model.

We would like to remark also some other points:

a) The model is conceived to describe two-particle Lagrangian dispersion. This does not have to be considered as a limitation since it can be verified (not shown) that also one-particle, or absolute dispersion, statistics is reproduced coherently with the theory of diffusion.

b) The statistical quantities characterizing the relative dispersion process are not altered if one considers the separation between particles that belong to different pairs. This assures that the properties of dispersion are common to a whole set of tracer particles and are not limited to the single pairs separately.

c) The kinematic eddies are long-living coherent structures, i.e. with slowly decaying Eulerian auto-correlations. It can be verified that letting the velocity amplitude of the modes evolve in time as a stochastic variable, e.g. via a time-correlated Langevin equation, does not change the Lagrangian characteristics of the dispersion process.

d) The kinematic model can be used either as independent model of an ideal turbulent flow or as sub-grid or, more in general, small-scale model for turbulent dispersion inside some more complex and realistic large-scale model of, e.g., atmospheric or oceanic circulation. At this regard, a two-dimensional version of the model, see Lacorata et al. (2014), can be adapted to improve the simulation of horizontal mesoscale dispersion over large domains.

e) The scaling properties of the model are not limited to the Richardson’s regime but, in principle, any arbitrary scaling law can be imported in the model, depending on the available experimental information about a given system. For modelling dispersion in the ocean upper layer, for example, the 2D version of the model can be set up according to the information coming from Lagrangian drifter data.

Acknowledgements.

GL would like to thank “RITMARE” (Ricerca ITaliana per il MARE) Reseach Program for providing financial support and scientific motivations for this work. AV thanks the Kavli Institute for Theoretical Physics, CAS, Beijing, for the kind hospitality at the Workshop ”Non-equilibrium processes at the nanoscales” (July-August, 2016) where this manuscript was completed.