Eulerian and Lagrangian studies in surface flow turbulence

Experimental and numerical studies of turbulent fluid motion in a free surface are presented. The flow is realized experimentally on the surface of a tank filled with water stirred well below the surface. Numerically, it is modelled by free-slip boundary conditions. The surface flow is unconventional: it is not incompressible and neither kinetic energy nor enstrophy is conserved in the limit of zero fluid viscosity and in the absence of external driving as is the case for incompressible two-dimensional turbulent flows. The dynamics of passive Lagrangian tracers that are advected in such flows are dominated by rapidly changing patches of the surface flow divergence. Owing to compressibility, particles form clusters within multifractal mass distributions. Also studied is the motion of pairs and triplets of particles. The mean square separation shows an extended range with a reduced scaling exponent in comparison with the classical Richardson value. Clustering is also manifest in strongly deformed triangles spanned within triplets of tracers.


Introduction
Few studies of turbulence have been more important than the paper by L F Richardson in 1926 on the rate of separation of particle pairs in a turbulent flow [1]. He argued that even in such flows, which are entirely different from that of uncorrelated Brownian particle motion, one could ascribe a turbulent diffusivity D turb to their relative motion. Since the velocity difference of particle pairs increases with their separation [2], D turb must likewise be an increasing function of their separation.
Richardson's arguments were supported by measurements of D turb ( ) that spanned five decades in separations . The data he compiled, including his own, could be fitted to the equation D turb ( ) ∝ (t) 4/3 . This equation, and the data supporting it, could fairly be said to mark the beginning of our present understanding of turbulence. Whereas Kolmogorov presented an analysis of turbulence [3] that underpins the present paradigm, that paper did not appear until 15 years later. Amongst the wealth of Kolmogorov's predictions is this same value of 4/3. In recent years, there has been a resurgence of interest in Richardson diffusion although laboratory experiments are few [4,5], because tracking particles is difficult without recently developed highspeed photography [4,6], high-speed particle detectors [7] or acoustic methods [8]. With respect to direct numerical simulations, even within the highest possible resolution for homogeneous isotropic turbulence, it is hard to obtain an extended range with Richardson scaling [9,10].
The earliest studies of turbulent diffusion were atmospheric, and the first ones were carried out by Richardson and his associates [1]. Inspired by Richardson's idea, Stommel used tracers to probe flow structures on the ocean surface 20 years later [11]. He simply seeded the surface with white paper sheets and took photographs of their motion from an aeroplane. Another 30 years later more sophisticated experiments, but in the same spirit (tracking floaters from satellites), have led to new insights on the mesoscale and large-scale structures in the upper ocean such as for the most famous example, the Gulf Stream [12]. Vertical cut through a velocity snaphot illustrating the mechanism for particle coagulation. Particles collect at downwellings (D) and spread away from upwellings (U). Only a part of the whole simulation domain is shown here.
Here we review and extend experimental and numerical studies in a very similar system, turbulence on the surface of a fluid volume. We will discuss properties of the flow within both the Eulerian and Lagrangian picture, respectively. Experimentally, this is achieved by small particles that float on a tank filled with turbulent water. The floaters sample the horizontal components of the velocity, v x (x, y, z = 0, t) and v y (x, y, z = 0, t). Since water is incompressible at all points (x, y, z), one has for the surface at z = 0 ∇ 2 · v = ∂v x (x, y, 0, t) ∂x + ∂v y (x, y, 0, t) ∂y = − ∂v z (x, y, 0, t) ∂z . (1) Hence the tracers on the surface probe a compressible two-dimenisonal system, even though all measured velocities are much less than the speed of sound. Figure 1 illustrates how the divergence arises. Sinks appear where the fluid is moving downward (D) and sources where it is moving upward (U). Particles, initially distributed uniformly, gather at the sinks and flee the upwellings. Whereas the erratically moving water molecules in the turbulent flow can come up to the surface and retreat from it, the floaters cannot. However, they can exchange both kinetic energy and squared vorticity with the water below them, hence dimensional arguments cannot be invoked to explain the velocity spectrum of the floaters or their relative diffusion. The flow differs from two-dimensional (2D) turbulence because neither kinetic energy nor enstrophy are inviscid invariants [13]. This removes any compelling reason to expect that the exponent of the relative dispersion, for example, should turn out to be close to 3, and indeed we find that it does not. The present system can be considered as a precursor of real compressible turbulence in which the pressure field would be an independent variable. It can be realized rather straightforwardly in the laboratory as well as in the simulations. Furthermore, the current system shows some similarities to inertia-generated clustering, as manifest in a rain storm where big clouds are formed. This recently became a subject of great practical and fundamental interest [14]- [16].
Although only turbulent flows are studied here, particles floating on a fluid in laminar motion would probe compressibility effects that can induce their coagulation as well. The convective rolls seen in (free-surface) Rayleigh-Benard convection would lead to a similar situation The particles which are 10 µm in diameter appear as white and were initially dispersed uniformly over the surface. This image was captured 100 ms later. The scale below the snapshot is in cm.
[17] but the dispersion would differ because the flow is spatially smooth which favours exponential separation [18]. The clustering effect, which will be discussed here, also leads to the formation of sharp fronts of particle density as seen in figure 2. As stressed in [19,20], the evolution of these fronts is illuminated by tracking the evolution of triplets of particles for the 2D case, as opposed to pairs alone. In the above-cited works for incompressible systems, it is found that triplets forming equilateral triangles, flatten with time. Here, as in the evolution of pairs, the experiments and the direct numerical simulations (DNS) are in qualitative agreement [21]. Nevertheless, there is not yet a way to obtain any of the results presented here from a theoretical model or from dimensional arguments.
The outline of the paper is as follows. In the next two sections we will discuss the experimental arrangement and the numerical method in brief. Afterwards we review the Eulerian properties of the flow, such as longitudinal structure functions, vorticity structure functions, the surface flow divergence and the compressibility ratio. In section 5 we turn to the Lagrangian studies. The particle distribution can be associated with a mass density, so investigations on its multifractality are presented. We report findings on the dispersion of tracer pairs, triplets and the particle density. We conclude with a brief summary and an outlook. Experimental apparatus for surface studies. Two mirrors and a cylindrical lens produce a sheet of laser light at the surface or in the interior. A camera positioned above records the images at a frame rate of ∼300 Hz.

Experimental arrangement
The measurements were made in a square tank (1 m × 1 m) filled with water to a depth of 30 cm. The set-up appears in figure 3 and a more in-depth description can be found in [22]. The beam from a diode-pumped solid-state laser is spread into a sheet, illuminating neutrally buoyant polystyrene particles in the turbulently stirred bath (stirring scheme not shown here) or alternatively hollow glass spheres floating on the surface. The fast overhead camera enables one to track the particles, with their position and local velocity recorded in a computer. The image acquisition rate was adjusted for image size, and was typically 300 frames s −1 .
Not shown in this figure is the equipment for generating the turbulence in the tank. A 8 hp pump recirculated water in the tank through pipes coming in from opposite sides of the tank in a plane located several cm above the bottom. These pipes delivered water to an array of 36 short vertical pipes capped with rotating nozzles that divert the flow horizontally. Each nozzle acted like a lawn water sprinkler, sending jets into the bulk. Strong turbulence is generated with this scheme, while at the same time keeping the wave amplitude small. This arrangement provides an unobstructed view of the water surface. Initial investigations were done with a vertically moving grid [6].
It was essential for the success of the experiment to keep the surface clean. Exposure to the air and outgassing of the plexiglass walls allows a film to form that interferes with the turbulence-driven motion of the floaters. Therefore, the surface was continuously vacuumed by a small auxiliary pump connected to a vertical pipe that comes up to the surface. Without this cleaning, tracers form a film that destroys the clustering apparent in figure 2. The floating particles used in the velocity measurements on the surface were hollow glass spheres with a radius a of 200 µm and highly buoyant particles of specific gravity 0.25 with a = 50 µm. The Stokes time, τ s = 2a 2 ρ p /(9νρ f ), was roughly 1 µs, so that the tracers easily follow the flow beneath them. Other particles, such as mushroom spores and talc, were also used successfully. The velocity measurements in the bulk were made using almost neutrally buoyant polystyrene spheres 10 µm in diameter.
This scheme, for generating turbulence with minimum wave production, enabled the experiments to reach a Taylor microscale Reynolds number Re of ∼150. With the rootmean-square (rms) velocity in the tank, the kinematic viscosity ν, and the Taylor microscale λ = (v 2 rms )/ (dv x /dx) 2 , this Reynolds number will be defined as The mean flow velocity in the tank is very close to zero, so the direction x is arbitrarily chosen. At Re = 140, the rms velocity fluctuations in the tank was approximately 3.5 cm s −1 . The Reynolds number as well as several other turbulent parameters are listed in table 1. To measure the strength and uniformity of the turbulence generated in the tank, the floaters were replaced by neutrally buoyant particles, enabling measurements of the velocity field by illuminating a plane below the surface with a sheet of laser light. The measurements in figure 2 were made by densely covering the surface with 10 µm floaters and tracking their dispersal in time; however this scheme turns out to be impractical for quantitative measurements of particle pair separation. The reason is that even small wave motion affects the illumination of surface particles, making it difficult to track the same particle for an extended period accurately. A preferred scheme which is used here for all the quantitative measurements is to map out the velocity by tracking as many particles as possible over a short interval of time. The velocity at every point and at each measuring time t is then determined by an average of the local vector field about the point (x, y); the local velocity measurements, v x (x, y, t), v y (x, y, t) are completed in a few ms.
Within the limits of the in-house software that analyses the data, one obtains a large number of particle 'tracks', on which our interest centres. The particles used for measuring the velocity field by a particle-tracking algorithm (which was written by M Rivera) were small enough to follow the local flow field. It is apparent that, even in a divergence-free flow, particles having appreciable inertia will also cluster [15]. In the computer simulations described here, the particles being tracked have zero inertia, and the inertia of the floaters in the experiments is likewise negligible. Figure 4 shows the spatial distribution of the floating particles at successive instants of time after they were uniformly distributed on the surface at t = 0. The particle positions in each image were deduced from measurements of the velocity field, as discussed above. We refer to this scheme as the tracking of 'virtual particles', and it is used for all of the Lagrangian measurements. The data in this figure span the interval 30 ms t 600 ms. Note the strong similarity between the true particle distribution measurements in figure 2 and the image at t = 300 ms.
In the Richardson diffusion measurements presented here, the mean square separation 2 (t) of the virtual particles was averaged over thousands of particle pairs in each run and the data were further averaged over 27 runs. For each run the trajectories were followed for up to 10 s at a data collection rate of 100 images s −1 .

Numerical simulations
The 3D Navier-Stokes equations for an incompressible Newtonian fluid are modelled numerically by a pseudospectral method using fast Fourier transformations and a 2/3 de-aliasing. The equations of motion are advanced in time by a fourth-order Runge-Kutta scheme. Periodic boundary conditions are taken in the lateral directions x and y. Free-slip boundary conditions are applied in the vertical direction z. They imply that the surface is perfectly flat in contrast with the experiments where vertical variations are always present although they are kept small here.
An appropriate volume forcing density which is added to the Navier-Stokes equations sustains the turbulence in the bulk of the turbulent flow. We used a volume forcing that injects energy at a fixed rate into the flow [23] or one that sustains a large-scale linear shear profile. More details can be found in [24]- [26]. The numerical grid resolution for most cases was N x × N y × N z = 256 × 256 × 65 points for a volume with an aspect ratio of L x : L y : d = 2π : 2π : 1. The spectral resolution criterion of k max η > 1 is satisfied. Up to N = 36 000 Lagrangian tracers are tracked simultaneously with the flow in the surface z = 0. Therefore the equations for i = 1, . . . , N must be solved. The surface velocity between the grid mesh is interpolated by bi-cubic splines [27]. Figure 5 shows a typical snapshot of the turbulent surface flow (red vector

Longitudinal structure functions
The floating particles allow one to reproduce the velocity field which is advecting them. Consequently, one can study some statistical properties of the Eulerian fields, such as moments of velocity increments. The closed circles in figure 6 show the second-order longitudinal velocity structure function which is defined as measured at the surface at a Taylor microscale Reynolds, Re = 93. The open squares correspond to measurements at a depth z = z 0 = 1 cm, where Re = 120. The brackets denote an average over the raw velocity field. Both sets of measurements correspond quite well to the scaling form D 2 (R) ∝ R ζ 2 , with a scaling exponent ζ 2 , which is close to the classical Kolmogorov value [2]. The computer simulations agree with the experiments at the surface which was shown in [6]. They enable one to look deeper into the behaviour of the nth-order moments of the velocity  [6,25]. The structure functions were evaluated always in a plane of constant z 0 . The turbulence was driven by a largescale shear in vertical direction. The statistical analysis was performed with about 100 independent snapshots of the turbulent field at the surface.
increments, D n (R), on the surface versus the bulk. The higher moments are defined similar to (4) and we had to apply extended self-similarity (ESS) due to the moderate Reynolds number [29]. Figure 7 shows how the moments n = 2 to 6 vary as a function of depth z = z 0 . Data up to the centre plane at z 0 /d = 1 2 are shown. One sees that sufficiently far from these surfaces, the moments fit well to the ESS values of ζ n known from homogeneous isotropic volume turbulence [2,29]. As the surface is approached, however, deviations increase, suggesting that the velocity field is more intermittent there according to the extended self-similarity analysis. This increase in the deviations δD n,3 (z) = D n,3 − n/3 from classical scaling of n/3 occurs at a value of z that is roughly ν/v rms , where v rms is the rms velocity in the bulk. D n,3 is the nth-order local ESS scaling exponent. In [25] we also investigated two flows at different Reynolds numbers and found such behaviour for both cases. Nevertheless, we expect that this larger magnitude of the deviations from the classical scaling in the surface can be influenced by finite Reynolds numbers or the particular kind of forcing.

Vorticity structure function
The second-order structure function of the vorticity field For these data the bulk turbulence was generated by an oscillating grid rather than the scheme discussed in section 2. Inset: numerical simulation results for S ω (R) at the free-slip surface, with Re 100.
which is defined as was also measured and calculated at the surface [6]. The vorticity data at z = 0 appear in figure 8. For the measurements, Re = 35, and for the simulations (inset), it is about 100. Other parameters appear in the figure legend. The turbulent outer scale L, defined as the integral of the velocity correlation function, is L = 1.7 cm. In the simulations S ω ∝ R 2 for R in the viscous range of turbulence. The vorticity structure function saturates roughly at L/2. This indicates that the spatial correlations C ω (R) = ω(x + R)ω(x) of the surface vorticity decay rapidly and hence we get saturation according to the relation S ω (R) = 2 ω 2 − 2C ω (R). The abscissa in the DNS is in dimensionless units R/η, where η = (ν 3 / ) 1/4 is the Kolmogorov dissipation scale. That small scale marks the crossover from rough to smooth flow. In the experiments, η was too small to measure directly, but ancillary observations suggest that η 0.1 mm [28]. It is apparent, then, that the two sets of measurements and the simulations span different ranges of R/η. Nevertheless, the qualitative behaviour is found to be the same for both investigations.
As we have discussed in [25], the non-vanishing surface flow divergence introduces an additional source for the vorticity component ω which prohibits the enstrophy from being an inviscid invariant. Therefore, it is also difficult to give dimensional arguments for a possible algebraic scaling relation in the inertial range. The run at Re = 145 was obtained with isotropic volume forcing. It was found that the PDF is slightly skewed [33].
Since water is incompressible, generating turbulence in the bulk inevitably produces waves on the surface for the experiment. These effects are absent in the DNS. The waves can similarly disperse surface particles, as is well known from previous works [30]- [32]. Experiments were conducted to establish that the amplitude of the surface waves is small enough to have negligible influence on the results presented here [6]. y, 0, t), which was found to fluctuate in time and space. Its probability density function (PDF) P ∇ 2 ·v is plotted in figure 9. The measurements were made at Re = 100 and 140 (left panel) [21]. The appreciable size of the seed particles made it impossible to directly measure velocity differences on a scale R less than the dissipative scale, η.

Divergence of the surface flow turbulence
The velocity divergence data were extracted from an ensemble of measurements made at many instants of time and many points (x, y) in the surface. Note that the mean value of ∇ 2 · v is close to zero and that this function has exponential tails out to almost four standard deviations, σ ∇ 2 ·v . Its measured standard deviation is about 10 s −1 at Re = 140. The width of the PDFs of the surface divergence is of the order of δv(η)/η, where δv(η) is the rms of the velocity difference at the scale η. Simulation results of P ∇ 2 ·v are shown on the right panel of figure 9 where the effect Table 2. Numerical simulation data which are taken at z = 0. Several kinds of forcing were used. Letter I stands for isotropic volume forcing that injects energy into the flow at a fixed rate (see e.g. [23]) and I + S for a combination with a volume forcing that sustains a mean streamwise flow with a constant shear in z (see e.g. [24]), respectively. Re is the Taylor microscale Reynolds number and C is the compressibility ratio (cf (7)). of different kinds of forcing are compared. Whereas the exponential shape of the PDF seems to be a robust feature, the slope in the tails seems to depend on specifics of the particular driving of the turbulence.

Compressibility ratio
The Eulerian compressibility is defined as [34] This expression can be applied at the surface, z = 0, and in any x-y plane of the bulk, where it will be non-zero likewise because of the 3D incompressible nature of the flow there. For the measurements and the simulations we obtain C 0.5 ± 0.1 on the surface. This result is very robust, being insensitive to changes of the way the turbulence is generated in the tank [6,33], to variations in the forcing spectrum for exciting the turbulence in the DNS, or to the aspect ratio (see table 2).
It is possible to re-express C by using free-slip boundary conditions as imposed in the simulations at z = 0. They are given as v z = 0 and ∂ z v x = ∂ z v y = 0.
If the flow is isotropic on the surface, the two-point correlation C i,j (R) ≡ v i v j , with the notation v i = v i (x) and v i = v i (x + R), cannot depend on any direction in the surface and thus [35] Since the velocity field on the surface is compressible, the divergence of the correlation tensor field leads to Due to 2D isotropy, the resulting vector can be expressed as v i ∂v z /∂z = d(R)R i /R, where d(R) is a scalar function. Thus one gets Using equation (11), one can show that the longitudinal and transverse second-order structure functions, exposed to 2D isotropy, can be written as where E = 1/2 v i v i is the 2D mean energy, i.e. i = j = 1, 2 and D T 2 (R) the transverse second-order structure function. We then take the limit R → 0 and for simplicity consider the direction of R along the y direction. Taking the direction of R along y, we expand equation (13) at the limit R → 0. Up to O(R 2 ), the left-hand side would be proportional to D T 2 (R) ∼ (∂v x /∂y) 2 R 2 . However, the expansion of the right-hand side of equation (13) By applying equation (14) one may rewrite equation (7) as Consequently, if the normalized cross-correlation function at the surface, i.e.
Thus, the small deviation of the compressibility from 1/2 is related to the small nonvanishing values of the cross-correlation. Of course, this is not the reason why the compressibility value should be 1/2 or even close to it, but at least it provides another consistency check by measuring the cross-correlations on the surface. In fact, the numerical value of the crosscorrelation tends to be very small. We have checked these values in a series of runs with different Taylor Reynolds numbers and two ways of forcing (see table 2). For higher Reynolds numbers all the instances give a value close to 1/2 and the cross-correlations yield small values at the surface. Further investigations on the effects of compressibility will be published elsewhere [37].

Richardson's pair diffusion theory
Turbulence will disperse a passive contaminant much faster than Brownian motion. If particles initially form a small cloud of mean square radius 2 (t) , the cloud will grow with ever increasing speed because the relative velocity δv(r) of particles of separation r(t) in the cloud, increases as r 1/3 , according to Kolmogorov [3]. It is this phenomenon that was addressed by Richardson in 1926 [1].
Richardson [1] suggested that the problem could be treated by a 1D diffusion equation, but replaced the diffusion constant with a turbulent diffusivity increasing with a power of r. In [1], he presented experimental observations that could be reasonably well fitted to the equation, D turb (x) = Dx α with α = 4/3. He inserted this result into the equation where ρ(x, t) can be regarded as a particle density or a PDF. The diffusion constant D is needed to make the equation dimensionally correct. Richardson solved this equation for the initial condition ρ(0, t) = δ(t) and found where K is a normalization constant. Evaluating the second moment, or variance of ρ(x, t), yields the result, with b = 3. This power-law variation should hold only in the inertial range for which x η [1]. Computer simulations in 2D and 3D incompressible systems [10,38] and recent experiments [4,39] are reasonably consistent with the above result, aside from a small intermittency correction which reduces the exponent slightly [10,40,41]. The experiments and simulations to be presented here, however, differ significantly from this, presumably because the floating-particle system is one with compressibility. by the lifetime τ η of eddies of the size η. The water tank measurements were at Re 100. The curve shows an average of an ensemble of 27 runs, each one tracking thousands of particle pairs having an initial particle separation, (0) = 3 mm. The largest separation was 20 mm, which still is within the integral scale L = 30 mm. Particle trapping and the finite area of the measurement region prevented following particle separations to even larger values. It should again be noted that the separations were deduced from measurements of the velocity field, not from the true tracking of the separation of particle pairs. Inertial effects, which also induce clustering [15], are virtually absent in the measurements and are not included in the simulations.

Pair dispersion at the surface
The experimental data are consistent with a change in slope from the value, b = 2 (which is the leading order in the initially exponential separation), to a number close to b = 1.65 [21,28]. To our knowledge, this is the first study of pair separation of compressible turbulence in a laboratory setting, although 2 (t) is frequently measured at the surface of the ocean, where the motion of the tracked floaters is also that of compressible flow [39].
The simulation data in figure 10 display an inflection point near t = 2τ η . Here the particles are initially separated by a distance of 1η. Lack of a wide inertial range prevents obtaining a quantitative measurement of this exponent. This difference between the DNS and the experiment remains an open point and we can only speculate about the reasons. The lower spatial resolution, small vertical extensions in the free surface (which are absent in the free-slip case) and the fact that the initial separation is an order of magnitude larger in the experiments [42] might cause the masking of such intermediate growth in the experimental data. The interesting feature of the data is the reduced value of b in the inertial range relative to the incompressible value. We suggest that the near-ballistic scaling in the inertial range results from the correlated motion of tracers within the clusters, where two processes are competing; the competition is between short-distance correlated motion within a particle cluster and decorrelation due to breaking of clusters when new fluid appears from below. The latter randomizing effect might be depleting the slope.
Measurements were also made on the dispersion of a patch of particles initially confined to a square 2 cm in size. Under the action of the turbulence, the patch of floating particles moved to cover the area of the tank, although their (virtual) positions were tracked only out to separations of the order of 10 cm. The data appearing in the left panel of figure 11 show the relative probability for particle separations at four different times following the initial (discrete and uniform) distribution at t = 0. Relevant experimental parameters appear in the figure caption. The distribution P( ) was obtained from an ensemble of 10 separate measurements temporally spaced by several seconds. For all members of the ensemble, the particles were initially distributed on a square array with a grid spacing of 0.45 mm.
After a short time (roughly 6τ η ), the distribution of particle separations is sharply peaked at a value close to 10 mm. This peak is of no fundamental interest; it merely arises because the patch of particles at t = 0 consists of equally spaced points spanning a square 2 cm on a side. The initial probability distribution must then be peaked at a value of lying between the closest particle spacing (0.45 mm) and the maximum spacing of √ 2 × 2 cm. At subsequent times the relative maximum of P( ) decreases as the distribution broadens. The most notable feature is the apparent existence of a plateau at large separations. The population shift to larger increases from t = 0.075 s (orange trace) to 0.6 s (red trace). This effect is apparent in a measurement made at t = 1.2 s (data not shown). By this late time, P( ) exhibits two highly populated regions, with the centre of the secondary population at equal to three-quarters of the entire 9 cm extension of the viewing region.
The right panel of figure 11 shows an expanded view of the P( ) at small . Here we see that the relative probability of particles being separated by very short distances increases with time. In fact, for the longest measured time (circles in red), the probability of zero separation is nearly 70% of the maximum value it reaches. In the interior of the system, the probability distribution would decrease monotonically for large separations. The relatively large number of particles separated by large distances might be indicative of spatially separated high-density clusters, which are apparent in figure 4.
Although Richardson's law is observed in 2D and 3D incompressible Navier-Stokes turbulence, his model, which is in essence a diffusion model at large times, relies on selfsimilarity and fast decorrelation of the flow. The simplifying assumptions clearly do not apply to the surface, as the flow exhibits finite Lagrangian correlation time and multi-scaling. Recently, theoretical studies of a Gaussian ensemble of rapidly decorrelated flows, so-called Kraichnan flows [43,44], have given new basic insights. The existence of long-time diffusive limit and of explosive separation are proven explicitly and the effects of compressibility, dimensionality and spatial roughness of the flow on dispersion were exploited analytically [34].
For a d-dimensional Kraichnan flow, for which D 2 (R) ∼ R 2ξ , the model predicts a critical compressibility given as C c = d/ξ 2 and separating two phases of particle dispersion. Within the so-called weak compressibility regime, C < C c , the particle trajectories separate explosively. In the strong compressibility regime, C > C c , the separation rate is diminished. The particle trajectories collapse and the compressibility dominates the dispersion [34,45].
The competition between the finite Lagrangian correlation time and the eddy turnover time has been addressed [46] in relation to the dispersion problem. In a nutshell, if the eddy turnover time asymptotically dominates the Lagrangian correlation time at large spatial separations, then the dispersion displays a long-time diffusive behaviour. In the opposite limit, i.e. when the Lagrangian correlation time is much larger than the eddy turnover time, the particle pairs feel a frozen flow at long-time asymptotics. The existence of two limits, the diffusive and the frozen one, is demonstrated in a synthetic flow [44], [46]- [48].
Whether those results are applicable for a time-correlated Navier-Stokes flow remains an open issue. It is not even clear that in the presence of time correlations the dispersion properties would enjoy as much generality as in decorrelated flows. Further measurements at higher resolutions are required to investigate the Lagrangian correlation time as well as its relative magnitude to Eulerian correlation time for varying degrees of compressibility. The limited inertial range in the numerical simulations does not yet permit such quantitative measurements.

Shape distortion of triangles spanned by particle triplets
A turbulent velocity field will distort the shapes and volumes of structures formed by three or more particles lying in a plane (in 3D, the minimum number of particles that will capture the effect of shear in all directions is four). Presumably, one is interested in the evolution of these shapes when the three points forming the structure are initially separated by lengths R that lie in the inertial range. In [21], we reported the evolution of tracer triplets that formed initially equilateral triangles. By neglecting the centre of mass, the relative position of the three points can be conveniently expressed by the following two vectors a 1 = (x 2 − x 1 )/ √ 2 and a 2 = (2x 3 − x 1 − x 2 )/ √ 6 [49]. The distortion of a triangle is measured by the ratio I 2 = g 2 2 /R 2 which relates its radius of gyration, R = a 2 1 + a 2 2 , to the smaller half-axis, g 2 , of the smallest ellipse that covers the tracer triangle. The ratio I 2 is defined as follows in 2D: Because the rate of separation of particle pairs increases with separation, a small fluctuation in the shear will subsequently drive three equally spaced particles, into an elongated structure, with the average elongation growing with time. As a result, one might expect I 2 to become small. In [19], it was argued that after a sufficiently long time, random small-scale fluctuations tend to drive I 2 towards the Gaussian equilibrium value, I 2G = (1 − π/4)/2 for the 2D case. In the 2D incompressible case, the experiments show that the Gaussian limit is achieved in roughly 20t/τ η [49]. It was asserted in [19] that coherent shear should dominate in the viscous subrange and give way to the randomizing effects of turbulent mixing at larger scales. In figure 12 we demonstrate this distortion by a sequence of snapshots for two particular tracer triples. One of the two triangles gets flattened completely while the second one is quasi-collinear.
The experimental data for I 2 are in the left upper box of figure 13 and the simulations are in the box to the right. We always started with equilateral triangles for which I 2 = 1/2, i.e. a 1 ⊥ a 2 and both sides have the same length.
The measurements and the simulations in figure 13 show a strong dip below the asymptotic value of I 2G and an extended minimum that exists for triangles initially having a radius of gyration significantly larger than the dissipative scale. This behaviour can be accounted for by strong clustering that is apparent in figure 2. The finite compressibilty of the flow causes a different temporal evolution of the shape measure. In incompressible turbulence, this effect is absent, and I 2 very soon attains its random Gaussian value I 2G . In the experiments, I 2 ultimately becomes slightly larger than the Gaussian value. The triangles in the computer simulations show an even stronger tendency to remain flat and appear to approach an equilibrium value smaller than I 2G . Self avoidance of the finite-size particles and possibly the particle tracking resolution used in the experiments force a limit on the minimum attainable value of I 2 and cause its later growth. In the numerical study (point) particles can become arbitrarily close and thus I 2 can remain at a small value well below Gaussian.
A dominant number of triangles remain distorted for very long times. This is illustrated in the two lower panels of figure 13. The joint PDF of the normalized vector magnitudes |a 1 |/R and |a 2 |/R is plotted there. After a short time (lower left panel), the triangles are only slightly distorted, which is reflected by weak scattering about the equilateral value, |a 1 |/R = |a 2 |/R. The right panel shows the triangles have become strongly distorted at a sufficiently later time and we observe two peaks (a 1 a 2 and vice versa).

Particle density
The number of Lagrangian tracers per area (which will be the grid mesh square for the following) defines a density field which is described by the continuity equation  is stretched and folded in a complicated fashion. Large deviation theory [50] can be applied to determine the asymptotic distribution of the contraction and stretching rates, which are nothing other than the so-called Lyapunov exponents. Quite generally, the entropy function that appears in the exponent of such PDFs is proven to be convex which causes low-order density moments to decay, but higher and negative orders to grow [51]. Recently, a relation between the statistics of stretching rates in a smooth flow and the dimension of the density field was derived [52]. Such general analytical results have yet to be extended to rough flows. The present surface flow contains both smooth and rough parts, respectively. Moreover, it is strictly non-Gaussian and time-correlated. All this makes the analysis of the advected density field very complicated as the flow is impermeable to theoretical models of advection in short-correlated Gaussian strains. In [26,33], we have discussed several methods to describe the temporal evolution of the formation of clusters, such as the particle density distribution, or the mass collected within clusters, exceeding a certain size. The outcome is that after an exponentially fast clustering at the beginning, a long transient takes over in which  still larger clusters are formed. Within the current experiments, we were not able to reach an eventual stationary state. Such slow convergences of density moments was also reported in [52].
To illustrate this feature we conducted an analysis of moments of mass density and the results are shown in figure 15. We observe a slight increase for the whole later stage of the R Figure 16. Scaling behaviour of the low-order mass density moments µ(B R (r)) q with respect to radius R. The values of q are −2, 0, and 2; 50 subsequent mass distributions are used for averaging and the time periods are indicated in the legend. The samples were separated by τ η /3 from each other. Inset: the spectrum GHP(q) for q = −2, . . . , 2 in steps of 0.015 is shown. For every particular q value we calculated the moment of mass density over R/η (as shown in the main figure) and evaluated the GHP spectrum by a fit over the range of the four smallest box lengths.
advection. Additionally, in the inset of figure 15 we plotted the evolution of global maximum of the particles per grid cell which does not saturate for t/T = 3. Resorting to the calculation of mass density moments , i.e. µ(B R (r)) q R ξ q and computing the mass µ(B R (r)) = ρ(r) dr on a ball B R (r) of radius R located at r we evaluated the Grassberger-Hentschel-Procaccia (GHP) spectrum [53,54]. It is given by Evidently the spectrum GHP(q) is related to the scaling exponent of the mass moment via q GHP(q + 1) = ξ q+1 . Moments of mass on a ball of radius R are calculated at the surface here (z = 0). As figure 16 shows, the non-linear character of ξ q as a function of q resembles qualitatively the intermittency of the mass fluctuations. However, the saturation level of the spectrum D q changes with time resembling the long transients again. This adds evidence that after an early fast stage of growth, the mass concentrates on smaller and smaller regions of the surface as additional accumulation continues slowly. A firm quantitative statement on the generalized dimensions would require far longer runs to achieve a stationary limit for the mass density moments which exceed our current numerical capabilities.

Summary
The relative motion of passive tracer particles that float on a turbulent body of fluid form a compressible system. Since they can freely exchange energy and vorticity with the incompressible turbulence beneath them, the floaters move very differently than in 3D or 2D incompressible systems. For one thing they seem to display stronger Eulerian intermittency than in the interior of the flow. For another, the turbulence creates convergent and divergent regions of surface particle density, and this affects their dispersion. While the floaters exhibit a significant range of spatial scales in which they separate algebraically, 2 ∝ t b , the exponent b is roughly half the Kolmogorov-Richardson value of 3. In addition, the surface compressibility alters the way that triangular structures on the surface are stretched by the underlying turbulence. These findings should be relevant to the motion of passive objects, such as buoys or pollutants moving on the sea [39]. Compressibility may also be important to take into account in the dispersion of active ocean surfactants, such as phytoplankton [55].