P Meunier

K Hourigan

Mixing in a vortex breakdown flow

This paper presents experimental and theoretical results on the mixing inside a cylinder with a rotating lid. The helical flow that is created by the rotation of the disk is well known to exhibit a vortex breakdown bubble over a finite range of Reynolds numbers. The mixing properties of the flow are analyzed quantitatively by measuring the exponential decay of the variance as a function of time. This homogenization time is extremely sensitive to the asymmetries of the flow, which are introduced by tilting the rotating or the stationary disk and accurately measured by Particle Image Velocimetry (PIV).

In the absence of vortex breakdown, the homogenization time is strongly decreased (by a factor 10) with only a moderate tilt angle of the rotating lid (of the order of 15 degrees). This phenomenon can be explained by the presence of small radial jets at the periphery which create a strong convective mixing. A simple model of exchange flow between the periphery and the bulk correctly predicts the scaling laws for the homogenization time.

In the presence of vortex breakdown, the scalar is trapped inside the vortex breakdown bubble, and thus increases substantially the time needed for homogenization. Curiously, the tilt of the rotating lid has a weak effect on the mixing, but a small tilt of the stationary disk (of the order of 2 degrees) strongly decreases (by a factor 10) the homogenization time. Even more surprising is that the homogenization time diverges when the size of the bubble vanishes. All these features are recovered by applying the Melnikov theory to calculate the volume of the lobes that exit the bubble. It is the first time that this technique has been applied to a 3D stationary flow with a non-axisymmetric perturbation and compared with experimental results, although it has been applied often to 2D flows with a periodic perturbation.

Introduction

It is well known that helical flows are unstable with respect to vortex breakdown, leading to a complex geometry of the streamlines with a barrier to transport between the inner part and the outer part of the bubble. The goal of this paper is to quantify experimentally the mixing in this 3D flow and to compare it with available theories derived for chaotic advection.

Vortex breakdown usually refers to a recirculating bubble that occurs past a stagnation point due to the rapid expansion of a thin vortex into a much broader one [START_REF] Leibovich | The structure of vortex breakdown[END_REF]. It is a surprising and practically important phenomenon, which has been observed in many different swirling flows from geophysical to industrial applications. Vortex breakdown was first observed over the delta wings of aircrafts [START_REF] Wentz | Vortex breakdown on slender sharp-edged wings[END_REF][START_REF] Hall | Vortex breakdown[END_REF][START_REF] Lowson | Vortex breakdown control by delta wing geometry[END_REF] where it creates a sudden drop of the lift and an increase in the drag, possibly leading to a loss of aircraft control [START_REF] Escudier | Vortex breakdown: Observations and explanations[END_REF]. By contrast, it is beneficial in the geophysical context since it decreases the maximal velocities of tornadoes and limits their destructive power [START_REF] Burggraf | Continuation or breakdown in tornado like vortices[END_REF][START_REF] Rotunno | The fluid dynamics of tornadoes[END_REF]. The structure and dynamics of these tornadoes is governed by the mixing of hot air within the core of the vortex. Artificial tornadoes within a chimney have been proposed as a way of converting green energy into electric power. But the poor understanding of the stability and mixing properties of these flows have prevented the development of an operational prototype. Vortex breakdown has also been used as a flame holder in combustion devices [START_REF] Gupta | Swirl flows[END_REF], where the mixing of fuel and oxygen is of primary importance. Finally, vortex breakdown flows are of interest for bioengineering applications since they can provide a smooth and efficient mixing without moving a propeller in the cavity [START_REF] Dusting | A fluid dynamics approach to bioreactor design for cell and tissue culture[END_REF][START_REF] Thouas | A bioreactor model of mouse tumor progression[END_REF]. This is beneficial for bioreactors where the growth of cells requires oxygen to be distributed as fast as possible, but where strong shears or moving objects should be avoided since they can damage cells. For such bioreactors, the goal is thus to maximize the mixing with a minimum velocity in the cavity, which will be the guiding line of the paper.

The structure and stability of a vortex breakdown flow has been actively discussed in the past century. Early experiments in a pipe [START_REF] Harvey | Some observations of the vortex breakdown phenomenon[END_REF][START_REF] Sarpkaya | On stationary and travelling vortex[END_REF][START_REF] Faler | An experimental map of the internal structure of a vortex breakdown[END_REF]) have revealed the presence of an axisymmetric bubble, which may also contain a single or a double helix (depending on the swirl and Reynolds number) with a strong hysteresis between these three regimes. The presence of the axisymmetric bubble has been explained theoretically by [START_REF] Benjamin | Theory of the vortex breakdown phenomenon[END_REF] and validated numerically by [START_REF] Wang | The dynamics of a swirling flow in a pipe and transition to axisymmetric vortex breakdown[END_REF] as a transition from a supercritical flow (without waves propagating upstream) to a subcritical flow (with waves propagating upstream) analogous to a hydraulic jump. [START_REF] Ruith | Three-dimensional vortex breakdown in swirling jets and wakes: direct numerical simulation[END_REF] further showed by Direct Numerical Simulations that helical disturbances constitute a secondary instability that grow over a finite time on the axisymmetric bubble. This secondary instability can be explained by a transition from a convective to an absolute instability [START_REF] Gallaire | Spiral vortex breakdown as a global mode[END_REF]) and leads to a global mode whose structure has been described by [START_REF] Meliga | A weakly nonlinear mechanism for mode selection in swirling jets[END_REF].

Vortex breakdown has also been observed in a simpler configuration consisting of a closed cylinder with one lid rotating [START_REF] Vogel | Experimentelle Ergebnisse über die laminare Strömung in einem zylindrischen Gehaüse mit darin rotierender Scheibe[END_REF][START_REF] Ronnenberg | Ein selbstjustierendes 3-komponenten-LDA nach dem Vergleighstrahlverfaren, angewendt für Untersuchungen in einer stationaren zylindersymmetrischen Drehströmung mit einem Rühckströmgebiet[END_REF]. This is the topic of this paper since the flow is very stable experimentally. [START_REF] Escudier | Observations of the flow produced in a cylindrical container by a rotating endwall[END_REF] showed experimentally that vortex breakdown appears in a finite band of Reynolds number if the height to radius ratio H/R is larger than 1.5. There is a single axisymmetric bubble until H/R = 1.95, where a second bubble can appear and even a third one around H/R = 3.35. In our experiment, the aspect ratio is equal to 2 such that the flow first bifurcates toward an axisymmetric time-dependent flow at a Reynolds number around 2700. Indeed, it was found numerically by [START_REF] Lopez | Axisymmetric vortex breakdown. part 3 onset of periodic flow and chaotic advection[END_REF] and theoretically by [START_REF] Gelfgat | Stability of confined swirling flow with and without vortex breakdown[END_REF] that periodic perturbations appear, leading to a vertical oscillation of the bubble. Note however that non-axisymmetric perturbations such as a double, triple or quadruple helix may appear before the periodic perturbations, but only for aspect ratios larger than 3.3 [START_REF] Gelfgat | Three-dimensional instability of axisymmetric flow in a rotating lid-cylinder enclosure[END_REF][START_REF] Lopez | Rotating and modulated rotating waves in transitions of an enclosed swirling flow[END_REF][START_REF] Sorensen | Multiple helical modes of vortex breakdown[END_REF][START_REF] Lopez | Three-dimensional swirling flows in a tall cylinder driven by a rotating endwall[END_REF].

The presence of a periodic perturbation superimposed on the axisymmetric vortex breakdown is interesting because it breaks the barrier to transport. Indeed, the axisymmetric streamlines surrounding the bubble are impermeable because they connect the two hyperbolic (or heteroclinic) points located on the axis. The transport of particles from inside to outside the bubble is thus completely governed by the additional periodic perturbation. [START_REF] Lopez | Axisymmetric vortex breakdown. part 3 onset of periodic flow and chaotic advection[END_REF] showed numerically that it advects the particles on the unstable manifold of the hyperbolic point. This is very similar to the theoretical result of [START_REF] Rom-Kedar | An analytical study of transport, mixing and chaos in an unsteady vortical flow[END_REF] on a 2D periodic flow of two oscillating point vortices.

They showed that the particles exiting the bubble are located in the lobes between the stable and the unstable manifold of the heteroclinic points, and that the area of these lobes can be calculated by the Melnikov function. The Melnikov theory is a classical tool of dynamical systems (see e.g. [START_REF] Guckenheimer | Non-Linear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF][START_REF] Wiggins | Introduction to Applied Nonlinear Dynamical Systems and Chaos[END_REF] which gives a simple formula for the "distance" between the stable and the unstable manifold.

However, the barrier to transport can also be broken by a symmetry breaking, as has always been found experimentally in the steady case. Indeed, there has been a long debate on the structure of vortex breakdown since it was different depending on whether it was predicted by axisymmetric numerical simulations or observed in the experiments that revealed asymmetric tongues of dye exiting from the bubble due to what was thought to be an intrinsic asymmetry of the bubble [START_REF] Spohn | Experiments on vortex breakdown in a confined flow generated by a rotating disc[END_REF]. It was finally proven [START_REF] Thompson | The sensitivity of steady vortex breakdown bubbles in confined cylinder flows to rotating lid misalignment[END_REF]) that these tongues can be due to very small imperfections of the experimental set-up, since a misalignment of 0.01 • of the rotation axis (which is of course inevitable) is sufficient to recover the large asymetries observed in the experimental images. This phenomenon highlights the high sensitivity of the flow to asymetries, which will be introduced on purpose in this paper in order to quantify their effect on the mixing. [START_REF] Sotiropoulos | Chaotic advection in threedimensional stationary vortex-breakdown bubbles: Sil'nikov's chaos and the devil's staircase[END_REF] showed numerically that it creates a rich dynamics such as KAM tori and stable islands inside the bubble, which have even been observed in experimental Poincaré sections [START_REF] Sotiropoulos | Experiments on lagrangian transin steady vortex-breakdown bubbles in a confined swirling flow[END_REF][START_REF] Fountain | Visualization of three-dimensional chaos[END_REF]. They also showed that the stretching and folding of the lobes creates a cantorus at the periphery of the bubble, which is responsible for a chaotic release of the particles initially located inside the bubble. This leads to a devil's staircase (with fractal steps) for the number of particles as a function of time. This has been confirmed by [START_REF] Brons | The influence of imperfections on the flow structure of steady vortex breakdown bubbles[END_REF] who also noted that the asymmetry can be reduced by introducing a small offset of the axis of the rotating lid to counter-balance its misalignment. However, there has never been a quantitative prediction of the escape rate from the bubble, as was done in [START_REF] Rom-Kedar | An analytical study of transport, mixing and chaos in an unsteady vortical flow[END_REF]. Moreover, all these studies only analyze the evolution of point tracers (position, spatial structure, fractal dimension...). They have always neglected the effect of diffusion which is able to break the barrier to transport and introduce a weak escape rate predicted theoretically by [START_REF] Brons | Dye visualization near a three-dimensional stagnation point: application to the vortex breakdown bubble[END_REF].

By contrast, many studies have been devoted to the advection-diffusion of a scalar in a homogeneous turbulent flow (see e.g. the review of [START_REF] Warhaft | Passive scalars in turbulent flows[END_REF]. The goal is to derive the spectra, the Probability Distribution Functions (PDF) and the variance of the scalar. In the inertial range the spectra scale as k -5/3 but they scale as k -1 at smaller scales (between the Kolmogorov and the Batchelor scale), as explained by [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in a turbulent fluid. Part 1. General discussion and the case of small conductivity[END_REF] and [START_REF] Kalda | Simple model of intermittent passive scalar turbulence[END_REF] in a stationary regime. The PDF reveal a high intermittency (non Gaussian), which can be explained by a simple model of stretched filaments [START_REF] Duplat | Mixing by random stirring in confined mixtures[END_REF][START_REF] Meunier | The diffusive strip method for scalar mixing in two dimensions[END_REF]. However, the most important result is probably the fact that the variance decays exponentially, as was shown numerically in a 2D almost time-periodic flow [START_REF] Pierrehumbert | Tacer structure in the large-eddy dominated regime[END_REF][START_REF] Antonsen | The role of chaotic orbits in the determination of power spectra of passive scalars[END_REF][START_REF] Thiffeault | A bound on mixing efficiency for the advection diffusion equation[END_REF][START_REF] Meunier | The diffusive strip method for scalar mixing in two dimensions[END_REF] and in a 3D flow [START_REF] Toussaint | A numerical Eulerian approach to mixing by chaotic advection[END_REF]. This can be explained in a bounded domain by the linearity of the advection-diffusion equation [START_REF] Haynes | What controls the decay of passive scalars in smooth flows?[END_REF]. However, [START_REF] Gouillart | Walls inhibit chaotic mixing[END_REF] showed experimentally that the walls can have a non-trivial effect and lead to "anomalous" power law rather than an exponential decay.

In this paper, we will study experimentally the mixing inside a highly heterogeneous flow. Section 2 is devoted to the experimental set-up. The flow is described in detail in section 3. The mixing properties of the flow are studied in the absence of vortex breakdown in section 4, and in the presence of vortex breakdown in section 5. Conclusions follow in section 6.

Materials and methods

Experimental setup

We wish to study the flow inside a circular cylinder with a rotating tilted disk. The set-up has been described in detail in [START_REF] Ismadi | Experimental control of vortex breakdown by density effects[END_REF] and will only be presented briefly here. Figure 1 shows the experimental set-up, consisting of a cylindrical Plexiglas container with internal radius R = 32.5mm, filled with water. This inner cylinder is placed inside an octagonal housing full of water with flat exterior faces to prevent the deformation of images. The flow under study is located inside the inner cylinder between two disks separated by a height H, with an aspect ratio H/R = 2 ± 1%. The top disk is rotated at an angular velocity Ω varying from 0.6 to 3.6 rad/s using a stepper motor and a worm wheel allowing a very smooth rotation of the disk at all speeds.

The goal of the paper is to study the mixing properties of the flow inside the inner cylinder. For this purpose, a dye mixture is injected at the centre of the bottom disk through a 0.5mm hole connected to a 1ml syringe. A syringe pump creates a continuous injection at a volumetric flow rate Q = 0.02 ml/min for 5 minutes. This duration of injection was chosen because it is small compared to the typical time of the experiments (of the order of a few hours) needed for the dye to spread within the whole cylinder. But this duration is long enough for the velocity of the injection jet (0.17cm/s) to be small compared to the velocity of the disk periphery varying from 2 to 12cm/s. The injected volume (0.1ml) is, of course, very small compared to the total volume of the cylinder (200ml). The mixture contains 0.033% of Fluorescein in weight diluted in water, together with 0.088% of alcohol in order to obtain a neutrally buoyant mixture. Indeed, [START_REF] Ismadi | Experimental control of vortex breakdown by density effects[END_REF] showed that a weak density difference of the order of 0.01% is sufficient to modify the threshold of the vortex breakdown.

The Fluorescein dye is illuminated by a blue laser so that it fluoresces very brightly with a green color. As illustrated in Figure 1(a), the laser is expanded into a vertical sheet that is carefully placed at the centre of the inner cylinder to visualise the distribution of dye in a longitudinal section. A digital camera records images similar to the one shown in Fig. 1(b). Although the blue component of the images contains strong reflections of the laser sheet, the green component only contains the image of the dye and its intensity is proportional to the density of the dye at this high dilution rate.

As will be shown later, the mixing properties of this flow strongly depend on the asymmetries of the experiment. We have thus tested two types of asymmetries: a tilt of the rotating disk and a tilt of the stationary disk. As shown in Fig. 1(b), the tilt of the rotating disk is obtained by cutting the disk at an angle α compared with the plane normal to the axis of the disk. This creates an oscillating motion of the periphery (in a given section) superimposed on the rotation of the disk. This does not correspond to the motion of a disk whose axis is tilted with respect to the axis of the cylinder. Six different disks with angles α equal to 0, 2, 5, 7.2, 9.8 and 14.6 degrees have been used, with an accuracy of less than 0.1 degree. The tilt of the stationary disk is obtained by placing a thin disk on a wedge at an angle β, which is varied from -2.6 to 2.6 degrees. The uncertainty on the angle β is relatively large even in the absence of a wedge.

The flow is also studied using Particle Image Velocimetry (PIV), by injecting small red fluorescent particles in the cylinder. A 4 mega-pixels PIV camera acquires images through a red filter in order to remove reflections of the laser. The images are treated using an inhouse cross-correlation algorithm [START_REF] Meunier | Analysis and optimization of the error caused by high velocity gradients in PIV[END_REF], which gives 60 by 60 vectors in the longitudinal section. The acquisition frequency is equal to 8 times the rotation frequency of the top disk in order to apply some phase averaging of the velocity fields.

Theoretical and numerical formulation

The flow under consideration depends on four main dimensionless parameters. The aspect ratio h = H/R is kept to a constant value of 2. The Reynolds number Re = ΩR 2 /ν, ν being the kinematic viscosity of the water, is varied between 700 and 3000. The angle α of the tilted rotating disk is varied between 0 and 15 degrees. The angle β of the stationary disk is varied between -2.6 and 2.6 degrees in the plane normal to the laser sheet. The Schmidt number Sc = ν/κ (κ being the diffusivity of the dye) is an additional parameter which is large in the experiment and equal to 2000.

The flow is governed by the Navier-Stokes equation together with the incompressibility condition, which are non-dimensionalised using the radius R and Ω -1 as length and time scales, leading to

∂u ∂t + u • ∇u = -∇p + 1 Re △u.
(2.1)

In cylindrical coordinates (r, θ, z), the boundary conditions are u = 0 at the periphery (r = 1) and at the bottom disk (z = βr cos(θ)) and u = re θ (where e θ is the orthoradial unit vector) on the top disk defined by z = h + αr cos(θ -t). These equations are easily solved numerically at these moderate Reynolds numbers using Comsol Multiphysics which is a finite element method using standard Lagrange elements P1-P1. The time-stepping is governed by an IDA solver (Implicit Differential-Algebraic) based on a BDF scheme (Backward Differentiation Formula) described by [START_REF] Hindmarsh | Suite of nonlinear and differential/algebraic equations solvers[END_REF]. At each time step, the system is solved using the direct linear solver for hollow matrices-Pardiso (www.pardiso-project.org). The flow is first calculated as a time dependent flow from t = 0 to t = 100, with u = 0 at t = 0. The solution at t = 100 is then used as an initial condition for a stationary solver. The mesh is refined in the boundary layer, leading to 14,882 elements in 2D and 511,000 elements in 3D.

The concentration of the dye is a passive scalar c advected by the flow according to the linear diffusion-advection equation

∂c ∂t + u • ∇c = 1 ReSc △c, (2.2)
This equation is not solved numerically because the Schmidt number is large such that dye structures are very thin. The calculation would need an increased spatial resolution at his large Peclet number (P e = ReSc ≃ 10 6 ), which would be very difficult in 3D. This is why the numerical simulations are only used to derive the characteristics of the flow and not to study the mixing problem itself.

Phenomenology of the flow

We first describe the helical flow inside the cylinder with a rotating disk, which is well known to exhibit vortex breakdown. The effect of the tilt of the rotating disk is analysed in detail to explain the mixing properties of this flow in the next sections.

Vortex breakdown in the presence of a tilted disk

The rotation of the top disk creates a thin rotating layer below the disk, which induces a radial motion due to the centrifugal force. This Von Karman pumping creates a secondary longitudinal recirculation inside the cylinder, which is very well visualised at moderate Reynolds number (see Fig. 2a) by the dye going up along the axis of the cylinder. This toroidal motion is very weakly affected by the tilt of the disk, which creates a weak periodic perturbation leading to a small undulation of the dye streakline below the top disk. The dye accumulates on the top disk and is then advected radially toward the corners by the centrifugal motion. At the corners, the dye predominantly goes down along the periphery of the cylinder due to the recirculation. But a small part is ejected radially inward due to small peripheral jets, which will be analysed in detail in section 3.3.

When the Reynolds number increases, the undulation of the streakline on the axis increases because the upward velocity decreases on the axis. At the critical Reynolds number where the velocity vanishes on the axis, the topology of the flow suddenly changes. This leads to the well-known vortex breakdown bubble, clearly visualised in Fig. 2(b).

Here, the aspect ratio h is equal to 2 such that the critical Reynolds number is equal to Re = 1430 ± 30, in excellent agreement with the literature (see [START_REF] Escudier | Observations of the flow produced in a cylindrical container by a rotating endwall[END_REF]. This critical Reynolds number is very weakly affected by the tilt of the disk: it seems to decrease by 20 due to the tilt of the disk, although this could be within the uncertainty.

Figure 2(c) shows that for a larger Reynolds number (Re = 1900), the vortex breakdown is larger and contains two bubbles. However, when the Reynolds number is increased further, the size of the bubble decreases and the bubble eventually disappears above Re = 3000 ± 50, as shown on Fig. 2(d). This is again in excellent agreement with the results of the literature obtained in the absence of tilt. It can be noted that the tilt of the disk creates a small undulation of the streakline just above the bottom disk, at the location where the bubble disappeared (i.e., where the vertical velocity is minimum).

Potential flow induced by the tilted disk

PIV measurements have been performed in order to quantify the meridional flow. Figure 3(a) shows the mean flow in a vertical section in the absence of vortex breakdown. The velocity field clearly shows the axisymmetric recirculation with a strong upward flow on the axis and a downward flow at the periphery. It can be noted that the radial flow is thinner and stronger at the top than at the bottom because it is contained within the layer. The periodic perturbation due to the tilt of the disk can be measured by calculating the difference between the phase averaged velocity and the mean velocity. This perturbation is plotted in Fig. 3(b) in the absence of vortex breakdown. It is clear that there is a strong flow from right to left, due to the tilted disk which pushes the fluid as it rotates. This flow is mainly potential because it is created by the motion of the disk normal to its surface. It can thus be calculated using the theory for potential flows, by assuming that the velocity u = ∇ϕ derives from the potential function ϕ. The potential function ϕ is searched as a sum of potential modes of amplitude a n :

ϕ = α n a n k n cosh(k n z) cos(θ -t)J 1 (k n r), (3.1)
where J 1 is the Bessel function of the first kind. Each mode is indeed a solution at any time because it satisfies ∇u = △ϕ = 0. The slip boundary condition (u r = 0) at the periphery (r = 1) imposes the choice of the axial wavenumber k n

J ′ 1 (k n ) = 0, (3.2)
such that there is a discrete set of wavenumbers, whose first terms are equal to k 1 = 1.8412, k 2 = 5.3314, k 3 = 11.7060 and k 4 = 14.863. The slip boundary condition at the tilted disk can be written:

u z = αr cos(θ -t) at z = h -αr sin(θ -t). (3.3)
This slip boundary condition is a priori valid only when the thickness δ of the boundary layer (which scales as Re -1/2 ) is much smaller than the vertical displacement of the tilted disk (which is equal to α at r = 1), i.e., when αRe 1/2 ≫ 1. However, it is also valid when αRe 1/2 ≪ 1 because in the frame of reference precessing with the disk (i.e., with the z axis normal to the disk), the Poincaré force is of order α and is thus negligible compared to the viscous terms governing the boundary layer which are of order Re 1/2 . In this frame of reference, the boundary layer is thus axisymmetric, meaning that the thickness δ is independent of θ, such that the slip boundary condition (3.3) is still valid but at a distance δ below the disk, i.e., at z = h -αr sin(θ -t) -δ.

Since δ and α are small in both cases, the boundary condition can be replaced at leading order by the same condition applied at z = h:

u z = αr cos(θ -t) at z = h, (3.4)
A solvability condition is then found by introducing the decomposition (3.1) into this boundary condition, multiplying by rJ 1 (k m r), and integrating r from 0 to 1 because 1 0 J 1 (k n r)J 1 (k m r)rdr = 0 if n and m are not equal. This leads to the value of the amplitude a n which can be simplified as The velocity is obtained by differentiating (3.1):

a n = 1 0 J 1 (k n r)r 2 dr 1 0 J 2 1 (k n r)rdr sinh(k n h) = 2 k n J 0 (k n )(k 2 n -1) sinh(k n h) . ( 3
u = n a n   J ′ 1 (k n r) cosh(k n z) cos(θ -t) -J 1 (k n r) cosh(k n z) sin(θ -t)/(k n r) J 1 (k n r) sinh(k n z) cos(θ -t)   . (3.6)
This velocity field is plotted in Fig. 4(a) and is similar to the perturbation field found experimentally and plotted in Fig. 3(b). Moreover, there is a good quantitative agreement of the amplitude of the flow, as shown in Fig. 4(b) where the radial velocity on the axis is plotted as a function of time. The experimental result clearly shows that the potential flow is independent of the Reynolds number. However, it seems that the experimental velocity is slightly delayed compared to the theoretical prediction. This may be due to nonlinear effects inside the boundary layer, which have not been taken into account here.

Peripheral jets induced by the tilted disk

As shown on Fig. 5(a), the perturbation flow (i.e. the difference between the phaseaveraged flow and the mean flow) contains almost no vorticity in the bulk of the cylinder, because the potential flow calculated in the previous section is irrotational. However, there is some vorticity created at the top disk which is advected to the periphery of the cylinder. It should be noted that this periodic vorticity perturbation is small compared with the mean distribution of vorticity (not shown here). The mean vorticity is antisymmetric whereas the periodic perturbation is symmetric. At the periphery, this symmetric vorticity corresponds to two helical vortices of opposite sign which are advected downward by the mean flow. These vortices are due to the oscillating motion of the disk at the corner, which creates some oscillating vorticity.

These counter-rotating vortices create small radial jets which are responsible for the tongues of dye which are visible in Fig. 2(a). The radial velocity of these jets is plotted as a function of z close to the cylinder in Fig. 5(b). It seems that the wavelength of the undulations (i.e. the pitch of the helical vortices) is independent of the Reynolds The strength of these jets decreases when the vortices are advected downwards due to viscous diffusion. Indeed, the attenuation is stronger at smaller Reynolds numbers: the oscillations disappear after 2 wavelengths at Re = 1000 whereas they are still visible after 4 wavelengths at Re = 2490. Despite the large uncertainty in the measurement, it is possible to fit the envelope of these oscillations by a decaying function of the form:

u envelope r = u max r e µ(z-h) (3.7)
The attenuation parameter µ is plotted in Fig. 6(a) as a function of the Reynolds number for 3 different tilt angles. It is quite clear that this parameter is fairly independent of the tilt angle α and decreases with the Reynolds number approximately as Re -1/2 . This is coherent with the fact that these jets extend on a larger height when the Reynolds number increases. The amplitude of these jets is characterized by u max r in (3.7) which is plotted in Fig. 6(b). It is surprising to see that the amplitude decreases with the Reynolds number. This is probably due to the fact that these jets come from an interaction of the viscous boundary layer (of thickness Re -1/2 ) with the tilt angle α. This would explain why the amplitude scales as αRe -1/2 . Using these two empirical laws, it is possible to give a general formula for the radial velocity of these jets at r = 0.9:

u r (θ, z, t) = 15αRe -1/2 exp[100Re -1/2 (z -h)] sin[18(z -h) + θ -t] (3.8)
where α is given in radians. This empirical formula will be used in the following.

To summarize, the presence of a tilt of the rotating disk does not modify drastically the mean flow, including the vortex breakdown bubble. However, it leads to a periodic potential flow in the bulk of the cylinder and creates periodic radial jets at the top corners which are advected downward by the mean flow.

Mixing in the absence of vortex breakdown

4.1. Mixing for a flat disk

Figure 7 shows the temporal evolution of the dye after injection on a long term. The dye is advected by the meridional recirculation on a torus and is thus homogeneous along the streamlines. The dye is mainly located on the axis, with a thickness which increases with time. This indicates that there is some diffusion across the streamlines, which allows to spread the dye over the whole volume and thus to decrease the intensity of the dye. However, this effect is weak because the dye is still inhomogeneous at t = 4000 (i.e. after 1 hour of experiment), as shown by Fig. 7(d) where the dye is located on the axis and at the walls. This leaves a large region without any dye, indicating that this flow is fairly inefficient for mixing.

In fact, there is a small asymmetry in the dye pattern (particularly visible on the axis of Fig. 7a), which is due to small imperfections in the experimental setup. This alternate chevron pattern is probably due to a non-zero angle β because this asymmetry is stationary with time. It is probably this asymmetry which creates a diffusion across the streamlines rather than the molecular diffusion whose characteristic time is equal to the Péclet number (P e = ReSc ∼ 10 6 ) in non-dimensional units. This already shows that the mixing properties of this flow are highly dependent on the possible asymmetries of the flow.

Mixing for a tilted disk

To study quantitatively the effect of the asymmetry, the rotating disk is tilted of an angle α = 5 • to surpass the effect of the imperfections of the set-up. Fig. 8 shows the temporal evolution of the dye at the same instants in time. The dye is still advected by the meridional recirculation, leading to a dense region of dye on the axis. However, the diffusion across streamlines is much faster, as the dye has spread over the whole volume at t = 1000 and is completely homogeneous at t = 4000.

This enhanced mixing is due to the tongues of dye which appear on the lateral walls due to the radial jets created at the top corners and advected along these walls (see section3.3). These tongues are advected by the flow to fill the bulk of the cylinder, leading to a rapid homogenization of the dye within the whole cylinder. This is clearly visible in 8(b).

Homogenization time

The mixing efficiency of the flow can be studied quantitatively by comparing the Probability Density Function (PDF) of dye concentration for two different tilt angles. These PDF are homothetic to the PDF of intensity P (I) (i.e. the normalized number of pixels with an intensity I) which are plotted in Fig. 9(a) because the dye concentration c is proportional to the image intensity I. For α = 0, the PDF contains a high and narrow peak at I = 10 corresponding to the regions without dye. But there is also a small and where • denotes the spatial average for r < 0.6 and z ∈ [0.1; 1.9]. The evolution of the standard deviation is plotted in Fig. 9(b) for two different tilt angles. It is clear that it decays exponentially in both cases after a small transient stage (although there is a saturation of the standard deviation at late stages for α = 5 • , due to the noise in the measurements). It is also clear that the decaying of the standard deviation is much faster for α = 5 • than for α = 0, which proves quantitatively that the mixing is strongly enhanced by the tilted disk.

The decay of the standard deviation is very well fitted by a law

σ 2 (t) = σ 2 0 e -2t/T + σ 2 noise (4.2)
which is plotted as lines on Fig. 9(b). This allows to measure accurately the homogenization time T (i.e. the decaying time of the standard deviation) for any experiment.

It should be noted that this homogenization time does not depend on the amount of dye injected (which only modifies the multiplying factor σ 0 ). It was also checked that it is independent of the region taken to make the spatial average at late stages. This can be understood by the fact that the dye concentration c is solution of a linear equation (2.2) which admits at any time an infinite set of eigen modes c n (r). For α = 0, the flow is stationary such that it is the same eigen mode at any time and the solution can be written:

c(r, t) = A n c n (r)e -t/Tn (4.3)
where -1/T n is the eigenvalue of the mode and A n is the amplitude of the mode at t = 0, which is given by the initial condition c(r, 0) = A n c n (r). It should be noted that T n is always positive because it is a diffusive process; it is actually equal to 2ReSc c n (r) 2 dr/ |∇c n (r)| 2 dr (see e.g. [START_REF] Toussaint | Spectral decay of a passive scalar in chaotic mixing[END_REF]. In an infinite medium, there is an infinite number of eigenvalues such that T n may not be bounded. This is the case, for example, in the absence of velocity where a blob of scalar converges to a Gaussian distribution whose variance decays as 1/t and not exponentially. However, in a confined geometry, the small-scale cutoff imposed by diffusivity reduces the problem to a finite-dimensional system, at least for 2D flows [START_REF] Haynes | What controls the decay of passive scalars in smooth flows?[END_REF]. This permits one to define the maximum decaying time T 1 , which corresponds to the dominant mode remaining at late stages since it is the least damped mode. This means that at late stages, σ is equal to the standard deviation of c 1 (r) over the region considered for the treatment of the images multiplied by A 1 e -t/T1 (if A 1 = 0). This demonstrates that the homogenization time T is independent of the region considered (unless c 1 (r) vanishes on this region). It also explains why the standard deviation σ decreases exponentially. For α = 0, the flow is not stationary but periodic. The same analysis can be done over each period, such that the result is still valid for integer values of time. The standard deviation is thus an exponential decay multiplied by a periodic function and -1/T 1 corresponds to the Floquet multiplier.

The homogenization time T is plotted in Fig. 10(a) as a function of the Reynolds number in the absence of vortex breakdown. It seems fairly constant around 1200 for a tilt angle α = 5 • . This can be understood by the fact that the radial jets are stronger but are damped quicker at smaller Reynolds numbers. These two effects annihilate each other which leads to the independence of T with respect to Re.

The homogenization time is plotted in Fig. 10(b) as a function of the tilt angle α for all Reynolds numbers (without vortex breakdown). It is clearly decreasing as a function of α, and seems to be well fitted to by a law scaling as α -1 . This shows that the homogenization time of this flow is strongly dependent on the tilt angle of the rotating disk. But it poses a problem for the axisymmetric configuration with α = 0. In fact, the homogenization time T saturates probably due to small imperfections in the symmetry of the set-up or at least due to the molecular diffusion which gives an upper bound equal to the Péclet number (Pe ∼ 10 6 at these Reynolds and Schmidt numbers).

A simple model

It is possible to give an estimate of the homogenization time by calculating the amount of concentrated dye which is ejected from the periphery to the bulk of the cylinder due to the radial jets. Using the empirical formula (3.8) of the radial velocity, the inward flux of fluid going from the periphery to the bulk can be calculated at r = 0.9 as

Φ = z=h z=0 max(-u r , 0)rdθdz ≃ 0.27α(1 -e -100Re -1/2 h ) (4.4)
This is valid for large Reynolds numbers when the radial jets are attenuated on a length larger than the wavelength of the oscillation. Furthermore, when h is larger than Re 1/2 /100 ≃ 0.3, the flux becomes independent of the Reynolds number, equal to 0.27 α. This explains why the homogenization time T is independent of the Reynolds number. This flux creates an exchange between the volume V p of the periphery at the concentration c p and the bulk with a volume V b at the concentration c b . It leads to two coupled equations for c p and c b , whose solutions are exponentially decaying toward the mean concentration with a characteristic time:

T = V b V p Φ(V b + V p ) ≃ 1.8h α (4.5)
Here, V b = 0.81πh and V p = 0.19πh correspond to the volumes inside and outside r = 0.9, and α is given in radians. This gives a correct scaling law for the homogenization time as α -1 but this prediction is 30 times below the experimental measurement plotted in Fig. 10(b).

It should be noted that the theory depends on the choice of the radius r at which the flux is evaluated. However, when r is close to the cylinder, the radial velocity vanishes such that Φ should be proportional to 1 -r. This seems to be the case experimentally for r ∈ [0.85; 1] despite strong uncertainties in the PIV measurements near the boundary. For r close to 1, the volume V p ≈ 2π(1 -r)h is proportional to 1 -r and much smaller than V b such that T ≈ V p /Φ is in fact independent of r. However, for r smaller than 0.85, the velocity is smaller than its linear variation in 1 -r. This would lead to an overestimation of the flux Φ in the theory which could explain why the homogenization time T is smaller than in the experiment.

Mixing in the presence of vortex breakdown

Homogenization time

Figure 11 shows the evolution of the dye in the presence of vortex breakdown. The images are completely different from the case without breakdown because the dye is trapped inside the bubble. This creates a sharp interface between the bubble with a strong concentration and the rest of the flow with a weak and homogeneous dye concentration. This clearly shows that the streamline which delimits the bubble is a barrier for transport since it prevents any advection across the streamline in the axisymmetric case.

As a consequence, the dye takes a very long time to escape the bubble, as shown in figure 11. Indeed, even at t = 4000 (Fig. 11d), the dye is still concentrated within the bubble although it would be completely homogeneous without vortex breakdown at this tilt angle (see Fig. 8d). It is curious to see that the dye is able to enter the bubble during the injection stage probably because the dye is extremely concentrated at the hyperbolic point. This allows for a sufficient flux of dye toward the bubble (due to the weak nonaxisymmetric perturbations) before the dye is spread in the bulk of the cylinder. The images are totally different if the dye is injected outside of the bubble (with a black bubble surrounded by a weak uniform dye) although it leads to the same homogenization time.

Figure 12 shows that the variance decays exponentially but with a characteristic time much larger than in the absence of vortex breakdown. This is clearly shown in Fig. 13(a) where the homogenization time T suddenly jumps by a factor 3 at the onset of vortex Moreover, the data are quite noisy, as if the mixing was governed by some imperfections of the experiment rather than the tilt angle α. Finally, it can be noted by comparing Figs. 11(a,b) that the asymmetry of the vortex breakdown bubble is stationary and that it does not rotate with the rotating lid. It is thus possible that the asymmetries of the bubble are due to a weak angle β of the stationary lid, which is well known to create large disturbances rather than the imposed tilt α of the rotating disk. In order to check this assumption, the tilt angle β has been varied by ±2.6 degrees with a constant tilt angle α of 5 degrees. Figure 14 shows the structure of the bubble for these two angles.

It is clear that the asymmetry of the flow is governed by the tilt angle β and not by the tilt angle α because the two images are symmetric to each other rather than identical. Furthermore, Fig. 15 shows a large variation of the homogenization time T by a factor 10 when the tilt angle β is varied by only 2 • . It seems to be correctly fitted to by

T = 1600 β -β 0 , (5.1)
where β is given in degrees and β 0 = -0.3 • is a small angle due to the imperfections of the experimental set-up. An accurate measurement of the shape of the cylinder (from the images) indeed shows that there is probably a small angle of the order of 0.2 ± 0.1 degree. It is surprising to see that a small tilt angle β of the stationary lid (of the order of 1 degree) has more influence on the homogenization time T than a larger tilt angle α of the rotating lid (of the order of 10 degrees). This result is in agreement with the numerical study of [START_REF] Brons | The influence of imperfections on the flow structure of steady vortex breakdown bubbles[END_REF] who showed that the residence time inside the bubble is extremely sensitive to the tilt angle β. They explained this sensitivity by a decorrelation between the unstable manifold of the lower hyperbolic point P and the stable manifold of the higher hyperbolic point Q. This creates some undulations of these manifolds, responsible for the tongues visible in the visualizations of Fig. 14, as described by [START_REF] Lopez | Axisymmetric vortex breakdown. part 3 onset of periodic flow and chaotic advection[END_REF]. In the following, we will try to quantify these tongues and the mixing they generate using the Melnikov theory in this 3D stationary flow, as was done by [START_REF] Rom-Kedar | An analytical study of transport, mixing and chaos in an unsteady vortical flow[END_REF] in a 2D flow with a periodic perturbation.

Melnikov theory on a model flow

We first need to model the flow as an axisymmetric bubble u 0 together with a perturbation flow βu 1 proportional to the tilt angle. The axisymmetric flow, presented in Fig. 16, exhibits a bubble located between two hyperbolic (and heteroclinic) points P and Q on the axis, separated by 2ρ ≃ 0.36 at this Reynolds number (ρ being the radius of the bubble). The angular velocity u θ /r varies between 0.1 and 0.2 between the lower part and the upper part of the bubble. However, this variation is smaller when the size of the bubble decreases because viscous effects tend to smooth the vorticity. As a consequence, the axisymmetric flow is modeled by the simplest normal form which presents two hyperbolic points P and Q around z 0 and has a constant angular velocity ω. It can be written in cylindrical coordinates as:

u 0 = -λr(z -z 0 ), ωr, λ((z -z 0 ) 2 -ρ 2 ) + ξr 2 .
(5.2)

In the following, z 0 will be taken equal to 0 which corresponds to changing the origin of the z-axis. Close to the critical Reynolds number, λ, ξ and ω are constant but ρ vanishes and it should be Taylor expanded to first order, leading to

ρ 2 = ζ(Re -Re c ).
The constants are measured using the numerical simulations at Re = 1900, which gives λ = 0.25, ξ = 0.17, ω = 0.15 and ρ = 0.18 (i.e., ζ = 7.2 × 10 -5 ). It is well known [START_REF] Brons | The influence of imperfections on the flow structure of steady vortex breakdown bubbles[END_REF]) that a small non-axisymmetric perturbation creates a difference between the unstable manifold w u (Q) (i.e. all the points coming from Q at t = -∞) and the stable manifold w s (P ) (i.e., all the points ending in P at t = +∞) In the frame of reference rotating with the bubble at ω, the trajectories remain in a meridional plane, such that the problem reduces to a 2D flow with a periodic perturbation (due to the non axisymmetric perturbation rotating at -ω). If these manifolds intersect in a point r i , they will intersect an infinite number of times, because a period 2π/ω later, r i is advected in r i+1 , which thus belongs to both manifolds w u (Q) and w s (P ) (see Fig. 16b). As shown in Fig. 16(b), this leads to an infinite number of intersection points which delimits lobes between the two manifolds. Moreover, the grey lobe located between r i and r i+1/2 and delimited by the two manifolds is advected a period 2π/ω later into the next grey lobe between r i+1 and r i+3/2 . Finally, as noted by [START_REF] Rom-Kedar | An analytical study of transport, mixing and chaos in an unsteady vortical flow[END_REF], the lobes come from the inside of the bubble (i.e., above w u (Q) at Q) at t = -∞ and are advected outside of the bubble (i.e. above w s (P ) at P ) at t = +∞. As a consequence, the volume of this lobe corresponds exactly to the quantity of dye that exits the bubble during one period. The goal of the calculation presented in the following is thus to calculate the volume of this lobe as a function of the non-axisymmetric perturbation.

The perturbation βu 1 due to the tilt angle has been evaluated numerically by subtracting the velocity u 0 obtained in a simulation at β = 0 from the velocity u obtained in a simulation at β = 2 • . However, since these calculations are three-dimensional and since the perturbation is small (because it is proportional to the small tilt angle β) there is significant noise in the data. The radial profiles of velocity are plotted in Fig. 17 at θ = ±45 • . The main effect of the tilt β is to create a non axisymmetric perturbation in the axial velocity w (solid line), related to a positive (resp. negative) jet at θ = -45 • (resp. θ = 135 • ). This can be modeled locally as a shear oriented at θ 1 , which leads to the velocity field for the perturbation:

u 1 = 0, 0, γr cos(θ -θ 1 ) , (5.3)
with θ 1 = 135 • and γ = 0.02 if β is given in degrees.

The Melnikov theory is applied in the appendix in order to calculate the distance between the stable and the unstable manifolds w u (Q) and w s (P ). For this purpose, the problem is solved in the frame of reference rotating with the fluid at ω such that the trajectories become two-dimensional. Moreover, we apply a change of variables from (r, z) into (s = πr 2 , z) so that the 2D flow becomes Hamiltonian. The volume of the lobes is then easily evaluated such that the homogenization time can be written analytically:

T = √ 2ρ 3 (λξ) 3/2 cosh(πω/2λρ) βγω(λ 2 ρ 2 + ω 2 ) .
(5.4)

The first consequence is that this homogenization time is indeed inversely proportional to the tilt angle β. This simply comes from the fact that the size of the lobes is proportional to β since the Melnikov theory is obtained to first order. The second consequence (less trivial) is that the homogenization time scales as T ∼ ρ 3 cosh(πω/2ρλ), i.e., a function of ρ, and it thus diverges when the size of the bubble vanishes. This is exactly what has been observed experimentally since T diverges when the Reynolds number gets close to the critical Reynolds number Re c . Of course, there is a saturation of T due to the molecular diffusion when the size of ρ is of the order of the diffusion length T /ReS , which leads to a diffusive time

T = ρ 2 ReSc ≃ 200(Re -Re c ).
This explains the rapid increase of the homogenization time at the onset of the vortex breakdown. Finally, it explains why the perturbation due to the tilt angle α of the rotating lid is less efficient than the angle β of the stationary lid. Indeed, the tilted rotating lid creates a non-axisymmetric perturbation which rotates at an angular frequency equal to 1 compared to the laboratory, i.e., at an angular velocity 1 -ω in the frame of reference rotating with the bubble (at ω). In formula (5.4), ω should thus be replaced by 1-ω when looking at the perturbation due to α. Since ω = 0.1, this creates a term cosh[π(1-ω)/2λρ] which is extremely large (10 12 times larger than for the tilt β). This explains why a large tilt α = 15 • of the rotating lid has almost no influence on the mixing across the bubble although a small tilt β = 1 • of the stationary lid has a large impact.

However, despite this qualitative agreement, the quantitative agreement is poor because the numerical values lead to T = 250/β at Re = 1900 (when β is in degrees), which is 6 times smaller than the experimental result. This could be due to the fact that the concentration within the exiting lobe is not exactly equal to the concentration within the bubble. Indeed, it contains some fluid coming from an incoming lobe at earlier time. This entanglement is known to create a fractal structure in the absence of molecular diffusion (see [START_REF] Rom-Kedar | An analytical study of transport, mixing and chaos in an unsteady vortical flow[END_REF][START_REF] Sotiropoulos | Chaotic advection in threedimensional stationary vortex-breakdown bubbles: Sil'nikov's chaos and the devil's staircase[END_REF]). This has not been observed here because the Schmidt number is not large enough. However, this effect leads to a weaker escape rate of the scalar which could explain the discrepancy between the experiment and the theory.

Conclusion

In this paper, we have presented an experimental and theoretical study on the mixing properties of a helical flow in a motionless cylinder with a rotating end. It was found that the axisymmetric flow is very inefficient for mixing such that any small asymmetry completely governs the homogenization of the dye inside the cylinder. Indeed, the homogenization time is always 2 orders of magnitude smaller than the theoretical diffusive time R 2 /κ, because small imperfections of the experimental set-up (here a small tilt of the bottom) create a convective mixing much larger than the diffusive mixing. In order to quantify the influence of the symmetry breaking on the mixing properties, the top ends have been tilted with controlled angles and the homogenization time has been measured accurately as the decaying time of the variance.

The mechanism for mixing strongly depends on the structure of the flow which is well known to exhibit a vortex breakdown bubble for a finite range of Reynolds numbers. In the absence of vortex breakdown, the tilt of the rotating disk induces some inward radial jets at the periphery of the cylinder, which are advected by the mean recirculation inside the bulk of the cylinder. As a consequence, the homogenization time is 3 orders of magnitude smaller than the axisymmetric diffusive time for a small tilt of only 10 degrees. The homogenization time is inversely proportional to the tilt angle and independent of the Reynolds number. These scalings can be explained by a simple model taking into account the flux of scalar between the periphery and the bulk of the cylinder.

In the presence of vortex breakdown, the scalar is trapped within the bubble such that the homogenization time is completely governed by the escape rate from the bubble. Curiously, the tilt of the rotating lid has a weak impact on the mixing. However, the tilt of the stationary lid has a very large influence. The homogenization time is inversely proportional to the tilt angle and seems to diverge when the size of the bubble vanishes. This surprising result can be predicted by using the Melnikov theory to calculate the volume of the lobes which exit from the bubble. As far as we know, this is the first comparison between an experiment and a theoretical derivation of a mixing rate using the Melnikov function.

It should be noted that this Melnikov theory has been usually applied for a 2D flow with a periodic perturbation. It is generalised here for a 3D axisymmetric stationary flow with a non axisymmetric perturbation. However, the angular velocity has to be assumed constant on the whole stable/unstable manifold in order to restrict the problem to a 2D flow in the rotating frame of reference. The result is thus only asymptotically valid for small bubble sizes (i.e., at the onset of vortex breakdown). The full solution could be calculated since the Melnikov theory has been extended to n dimensions by [START_REF] Gruendler | The existence of homoclinic orbits and the method of melnikov for systems in rn[END_REF] and applied to a vortex breakdown flow by Holmes (1984). However, their calculation is more difficult because the homogeneous problem ξ = Du 0 (t)ξ has to be solved at any time.

Moreover, although this type of dynamics of a perturbation added to heteroclinic points has been extensively studied in the absence of diffusion, there has never been any result characterizing the effect of the diffusion. However, this effect might largely influence the mixing properties. For example, it was shown here that the fractal properties of the entangled lobes is not recovered because the diffusion is too large. It would be interesting to quantify the cut-off at small scale due to diffusion as a function of the Schmidt number. Diffusion might also smooth the devil's staircase which is obtained in the number of particles trapped in the bubble. It might even change the escape rate of the scalar. Indeed, for small diffusion, the time for the scalar to exit the bubble is probably governed by the diffusion outside of a KAM torus rather than outside of the bubble since the separatrix is indeed permeable due to the perturbation. In fact, it is clear that there is a gap between the theories developed for homogeneous turbulence with molecular diffusion but without flow structures and the theories developed for dynamical systems with inhomogeneous flow structures but with no diffusion. This paper is an attempt to reconcile the two domains. However, new models still need to be developed in order to predict properly the mixing properties of flows with heterogeneous structures. trajectory (s 0 (t -t c ), z 0 (t -t c )) to first order and leads to the simple result that it is equal to:

M (t c ) = β +∞ -∞
u 0 [s 0 (t -t c ), z 0 (t -t c )] × u 1 [s 0 (t -t c ), z 0 (t -t c ), t]dt.

(6.4)

where the cross product u 0 × u 1 is equal to u 0 s u 1 z -u 0 z u 1 s . In our case, this formula is simply written (using a change of variable τ = t -t c ): since dl = u 0 dt. Here t i , t i+1/2 are two consecutive times at which the manifolds intersect (i.e., where M (t 0 ) = 0); for example t i = 0 and t i+1/2 = π/ω. It should be noted that this result was derived in 2D by [START_REF] Rom-Kedar | An analytical study of transport, mixing and chaos in an unsteady vortical flow[END_REF]; however, it can be used in 3D here because the volume element 2πrdrdz is equal to the surface element dsdz by definition of s. The volume of each lobe is thus equal to:

M (t c ) = -β τ =+∞ τ =-∞
V L = 4 √ 2βγπ 2 (λρ 2 + ω 2 ) 3(λξ) 3/2 cosh(πω/2λρ) . (6.8)

Finally, after each rotation period 2π/ω, each lobe has been translated from r i to r i+1 , such that the rate of scalar c which exits the bubble is equal to cV L ω/2π. Since the volume of the bubble is approximately equal to V B = 4πρ 3 /3, the conservation of the scalar leads to V B dc/dt = V L ω/2π such that the scalar c decreases exponentially with a decaying time equal to:

T = 2πV B ωV L = √
2ρ 3 (λξ) 3/2 cosh(πω/2λρ) βγω(λ 2 ρ 2 + ω 2 ) . (6.9)

Figure 1 .

 1 Figure 1. Schematic of the experimental set up.

Figure 2 .

 2 Figure 2. Dye visualisation of the flow inside the cavity at the end of the injection phase. (a) Re = 1000, (b) Re = 1470, (c) Re = 1900 and (d) Re = 3500. α = 5 • , h = 2.

Figure 3 .

 3 Figure 3. (a) Mean flow and (b) periodic perturbation flow obtained by PIV in a longitudinal section perpendicular to the plane of the tilt angle. Re = 1000, α = 5 • , h = 2

Figure 4 .

 4 Figure 4. (a) Periodic potential flow obtained theoretically using (3.6). (b) Amplitude of the periodic flow on the axis at z = 1.8 obtained experimentally at Re = 1000 (•), Re = 1255 (•), Re = 1490 (×), Re = 1750 ( ), Re = 2000 (⋆), Re = 2250 (+), Re = 2490 (♦) and Re = 3036 (▽). α = 5 • and h = 2.

Figure 5 .

 5 Figure 5. (a) Periodic perturbation vorticity obtained by PIV at Re = 2000. (b) Vertical profile of radial velocity at r = 0.9 obtained experimentally at Re = 1000 (solid line), Re = 1490 (dashed line), Re = 2000 (dash-dotted line) and Re = 2490 (dotted line). α = 5 • , h = 2.

Figure 6 .

 6 Figure 6. (a) Attenuation parameter µ and (b) amplitude umax used to fit the envelope of the radial velocity at r = 0.9 using (3.7). Tilt angles are equal to α = 5 • (▽), α = 10 • ( ) and α = 15 • (•).

Figure 7 .

 7 Figure 7. Dye visualisation at t = 0, 300, 1000, 4000. Re = 1000, α = 0 • . (negative image of the green component)

Figure 8 .Figure 9 .

 89 Figure 8. Dye visualisation at t = 0, 300, 1000, 4000. Re = 1000, α = 5 • . (negative image of the green component)

Figure 10 .

 10 Figure 10. Decaying time as a function of (a) the Reynolds number at α = 5 • and (b) as a function of the tilt angle α for all Reynolds number without vortex breakdown.

Figure 11 .

 11 Figure 11. Dye visualisation at (a) t = 0, (b) t = 300, (c) t = 1000 and (d) t = 4000. Re = 1900, α = 5 • .(negative image of the green component)

Figure 12 .Figure 13 .

 1213 Figure 12. Temporal variation of the variance of the image intensity σ in the presence of vortex breakdown (Re = 1900). The solid line corresponds to the fit (4.2). α = 5 • .

Figure 14 .

 14 Figure 14. Dye visualisation at t = 0 for (a) β = 2.6 • and (b) β = -2.6 • . Re = 1900, α = 5 • .(negative image of the green component)

Figure 15 .Figure 16 .

 1516 Figure 15. Decaying time as a function of the tilt angle β of the bottom stationary disk for Re = 1900. α = 5 • .

Figure 17 .

 17 Figure 17. Numerical profile of the perturbation velocity u 1 = uu 0 as a function of r along a line oriented at θ = -45 • (a) and θ = +45 • (b) at z = -0.5 and β = 2 • . Dashed lines correspond to ur, dash-dotted lines to u θ and solid lines to uz.Re = 1900

  2λs 0 (τ )z 0 (τ )γ s 0 (τ ) π cos[ω(τ + t c )]dτ. (6.5)This integral can be evaluated by splitting the cosine function and using the residue theorem, which leads toM (t c ) = 2 √ 2βγπ 2 ω(λρ 2 + ω 2 ) sin(ωt c ) 3(λξ) 3/2 cosh(πω/2ρλ) . (6.6)As noted by[START_REF] Rom-Kedar | An analytical study of transport, mixing and chaos in an unsteady vortical flow[END_REF], this function completely governs the transport across the streamline because the volume of the lobes is equal to

Acknowledgements

The authors would like to thank Andreas Fouras for providing the original experimental set-up, some acquisition material and a proper space in his laboratory.

The authors also acknowledge financial support from Australian Research Council (ARC) grants DP0877327 and DP110102141.

Appendix

The purpose is to evaluate the volume V l of the lobes between the unstable and the stable manifolds of the hyperbolic points P and Q in the flow u 0 + βu 1 given by (5.2) and (5.3). First, we need to study the problem in the frame of reference rotating at ω such that the trajectories remain in a meridional plane, leading to a 2D problem with a periodic perturbation u 1 z = γr cos(ωt). Second, we make a change of variable because the Melnikov theory can only be applied for a flow which is Hamiltonian to first order. We thus use the surface s = πr 2 instead of the radial coordinate. Using the relation ds = 2πrdr, the trajectories (s(t), z(t)) are solution of the flow:

To first order, this flow is Hamiltonian because

The Hamiltonian (which is called the streamfunction in fluid mechanics) defined by u s = -∂H/∂z and u z = ∂H/∂s is equal to

The Hamiltonian is equal to zero on the streamline joining the two hyperbolic points P and Q, such that s = -2πλ(z 2 -ρ 2 )/ξ. Introducing this relation in the flow to first order leads to the differential equation ż = λ(ρ 2 -z 2 ), which is easily integrated, giving the trajectories to first order:

3)

The Melnikov function M (t c ) can now be calculated and corresponds to the distance D(t c ) between the stable manifold of P and the unstable manifold of Q, multiplied by the norm of the velocity (u 0 s , u 0 z ). This function is evaluated at t = 0 and at the point (s 0 (t c ), z 0 (t c )) where t c can be varied from -∞ to +∞. The Melnikov theory (see [START_REF] Guckenheimer | Non-Linear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF], for a demonstration) can be applied by following the


Alain Bretto

Alain Faisant

Luc

come

New graphs related to (p,6) and (p,8)-cages


Alain Bretto

published or not. The documents may come


Pierre Tandeo
email: pierre.tandeo@telecom-bretagne.eu

Ronan Fablet

René Garello

Spatio-temporal segmentation of mesoscale ocean surface dynamics using satellite data

Multi-satellite measurements of altimeter-derived Sea Surface Height (SSH) and Sea Surface Temperature (SST) provide a wealth of information about ocean circulation, especially mesoscale ocean dynamics which may involve strong spatiotemporal relationships between SSH and SST fields. Within an observation-driven framework, we investigate the extent to which mesoscale ocean dynamics may be decomposed into a superimposition of dynamical modes, characterized by different local regressions between SSH and SST fields. Formally, we develop a novel latent class regression model to identify dynamical modes from joint SSH and SST observation series. Applied to the highly dynamical Agulhas region off South Africa, we demonstrate and discuss the geophysical relevance of the proposed mixture model to achieve a spatio-temporal segmentation of the upper ocean dynamics.

I. INTRODUCTION

In the last two decades, multi-satellite measurements of altimeter-derived Sea Surface Height (SSH) and multi-sensor measurements of Sea Surface Temperature (SST) have provided a wealth of information about ocean circulation and atmosphere-ocean interactions. As a depth-integrated quantity dependent upon the density structure of the water column, altimeter SSH estimations capture mesoscale structures, horizontal scales of 50 km to few hundred kilometers, and allow to retrieve surface currents using the geostrophy balance. This emerging and rich mesoscale circulation further stirs the largescale SST fields. Accordingly, our picture of upper ocean dynamics has considerably evolved towards a complex system characterized by strong interactions, whose spatio-temporal variability extends over a wide range of scales. Furthermore, several studies (see e.g. [START_REF] Lapeyre | Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory[END_REF] or [START_REF] Isern-Fontanet | Three-dimensional reconstruction of oceanic mesoscale currents from surface information[END_REF]) rationalize and demonstrate that fields of SST can become an active tracer coupled to the dynamics leading to strong correlations with SSH fields. To these assumptions the upper ocean dynamics may be simply predicted from surface density horizontal variations possibly dominated by SST variations. For such a case, a linear transfer function shall be identified between SSH and SST fields (cf. [START_REF] Isern-Fontanet | Potential use of microwave sea surface temperatures for the estimation of ocean currents[END_REF]).

Within an observation-driven framework, one may consider joint PCA/EOF (Principal Component Analysis, Empirical Orthogonal Functions) procedures to decompose the relationships between SST and SSH fields, as widely used in ocean sensing applications (cf. [START_REF] Leuliette | Coupled pattern analysis of sea surface temperature and TOPEX/Poseidon sea surface height[END_REF], [START_REF] Casey | Sea surface temperature and sea surface height variability in the North Pacific Ocean from 1993 to 1999[END_REF]). Such EOF-based schemes would however resort to a single linear and global model. As such, this model could not address the spatial non-stationarity of the SST-SSH relationships. By contrast, we here consider local linear transfer functions. We assume that locally, upper ocean dynamics may be analyzed according to a finite mixture model, where each component of the mixture is characterized by a local SST-SSH linear transfer function. This mixturebased representation relates to a nonlinear model. In this paper, we propose to investigate such a model to (i) develop a probabilistic learning-based setting for the inference of such mixture models and the spatio-temporal segmentation of the identified dynamical modes (i.e., the different components of the mixture model), and (ii) evaluate the extent to which such mixture models are geophysically relevant to characterize the upper ocean dynamics over active ocean regions.

Hereafter we consider the Agulhas region. The paper is organized as follows. Section II presents the remote sensing data. Section III describes our probabilistic learning-based model. In Section V, the application to satellite observations is evaluated. We further discuss and summarize the key results of our investigations in Section VI.

II. DATA

As SSH and geostrophic surface current (U,V) data, we use the daily delayed time Maps of Absolute Dynamic Topography (MADT) produced by Collecte Localisation Satellites (CLS) available online at http://www.aviso.oceanobs.com/. This information combines the signal of several altimeters onto a 1/3 degree Mercator projection grid. We use the 2004 data since four altimeters were available (Jason-1, Envisat or ERS-2, Topex/Poseidon and GFO). As SST data, we use optimally interpolated microwave SSTs provided by Remote Sensing System (RSS) available online at http://www.ssmi.com/. It combines the signal of three microwave radiometers (TMI, AMSR-E and WindSAT) which are robust to the presence of clouds. The spatial resolution is 1/4 × 1/4 degrees and the temporal resolution is the same as the MADT data, i.e. daily. We bilinearly interpolate the MADT data onto the SST grid. We focus on the Agulhas region between longitudes 5 • E to 65 • E and latitudes 30 • S to 48 degreeS. An example of SST and SSH fields with the associated geostrophic currents are given in Fig.

III. METHODOLOGY

A. Patch-based approach

As mentioned above, recent theoretical and numerical experiments have stressed that upper ocean dynamics may be characterized by couplings between SSH and SST according to the following relationship in the Fourier domain (cf. [START_REF] Held | Surface quasi-geostrophic dynamics[END_REF]):

F H ( SSH) = -γ|k| -α F T ( SST) ( 1
)
where k is the horizontal wavenumber vector, F H and F T are linear filters of SSH and SST respectively. The α parameter sets up the effective coupling between surface fields. For α = 1, Eq. ( 1) resorts to the Surface Quasi-Geostrophy (SQG) model. In [START_REF] Isern-Fontanet | Potential use of microwave sea surface temperatures for the estimation of ocean currents[END_REF], F H and F T were were band-pass filters between 80 km and 300 km. As α increases, the smoothing increases and couplings decreases. For α = 2, the filtered SST would trace the vorticity. Formally, Eq. ( 1) states that surface currents can be regarded as spatial derivatives of a filtered version of the SST field. The parameter γ relates to a normalization constraint. In general, parameters γ and α as well as the definition of the filters F H and F T , may spatially vary such that a single linear transfer function as in Eq. ( 1) is unlikely to apply globally. These considerations led us to hypothesize that zonal and meridional geostrophic surface currents (U,V) and SSH can still locally relate to SST derivatives, but according to a finite set of K linear transfer functions, hereafter referred to as K dynamical modes. Formally this is stated in the Fourier domain as

SSH = H k ( SST) (2)
where H k now characterizes the k th dynamical mode, which locally relates SSH and SST fields through linear filter H k .

Using a matrix formulation, Eq. ( 2) is rewritten in the real domain as a patch-based linear regression

Y(s, t) = H k (s, t)X(s, t) (3)

B. Latent class regression model

Our objective is to identify K different hidden surface dynamical modes from a joint set of SST patches (p-dimensional vector X) and SSH with the associated zonal and meridional surface currents (3-dimensional vector Y) as illustrated in Fig. 2. We assume that the conditional likelihood of Y given X resorts to a mixture of Normal distributions such as

p (Y|X, θ) = K k=1 λ k N k (Y; Xβ k , Σ k ) (4)
where N k represents a multivariate Gaussian probability density function with mean Xβ k and covariance Σ k , and λ k is the prior probability of mode k. In the literature, this model is referred to as a "latent class regression" or "clusterwise regression" (cf. [START_REF] Desarbo | A maximum likelihood methodology for clusterwise linear regression[END_REF]). To simplify the notations, we store the overall parameters of the model in θ = (λ 1 , . . . , λ K , β 1 , . . . , β K , Σ 1 , . . . , Σ K ) and estimate them using a classical maximum likelihood criterion (see [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] for more details). Then, we exploit the inferred model with parameters θ to perform a spatio-temporal segmentation of the underlying dynamical modes. More precisely, for any spatial location s and time t, using the Bayes' theorem, we evaluate the posterior likelihood that the dynamical mode is of type k such as

πk (s, t) = λk N k Y(s, t); X(s, t) βk , Σk p Y(s, t)|X(s, t), θ , ∀k. (5)
IV. RESULTS

A. Spatio-temporal segmentation

From the posterior likelihoods πk , we determine the segmentation maps of each dynamical mode as illustrated in Fig. 3 for a given date. The animations of the time series of these daily maps in the Agulhas current over 2004 are available online at: http://tandeo.wordpress.com/communications/articles/. A qualitative analysis of these maps of posterior likelihoods highlights a clear spatio-temporal segmentation of the different dynamical modes that can be interpreted from a geophysical perspective in terms of different geophysical processes.

V. CHARACTERIZATION OF OCEAN SURFACE DYNAMICS

For each dynamical mode, we report the observed distributions of current, height and temperature values (cf. Fig. 4). The first dynamical mode (red) characterizes very strong current magnitude and warm waters. It is primarily associated with the main Agulhas current that flows down the East coast of Africa through the Agulhas ridge. This mode also involves mesoscale eddies, the so-called warm core Agulhas rings with strong surface currents, low temperature gradients and middlerange SSH values around 0.5 m. The later seems to be a discriminative feature of this first mode. The second dynamical mode (green) mainly relates to the eastward Agulhas return current that hits a part of the South Atlantic current. It creates a subtropical front varying from 36 • S to 44 • S with strong eastward currents, middle-range SST gradients and large SSH values (about 1 m). The third (cyan) and fourth (blue) dynamical modes correspond to weaker surface currents. The third one is characterized by mid-temperatures and westward currents whereas the fourth one involves colder temperatures and eastward currents. Let us stress that the third dynamical mode involves large SST gradients but weak surface currents. In this mode, the SST can be clearly identified as a passive tracer of the surface upper ocean dynamics.

VI. CONCLUSION AND PERSPECTIVES

In this paper, we proposed an observation-driven framework to identify, characterize and track ocean surface dynamical modes. We rely on a latent class regression model, where the dynamical modes are characterized by a local linear Fig. 3. Maps of the posterior likelihoods given in Eq. ( 5) of the dynamical modes given the SST and SSH fields, the 1 st of September, 2004 within the Agulhas current. We use a four-class latent regression model fitted from the whole year 2004 (see text for details). For a given location and time, the sum of the of the four probabilities is equal to 1. In the subsequent, colors red, green, cyan and blue respectively distinguish the first, second, third and fourth dynamical modes. transfer function between SST, SSH and surface geostrophic current (U,V). This probabilistic approach locally relates the distribution of the SSH and sea surface currents conditionally to the SST via a nonlinear model: a Gaussian mixture of linear transfer functions. The statistical parameters of the model are estimated using a maximum likelihood approach.

We applied the proposed methodology to the 2004 daily 1/4 • × 1/4 • satellite SST and SSH image series. The reported results retrieved a relevant spatio-temporal decomposition of ocean surface dynamics in the Agulhas region according to four dynamical modes: (i) the main Agulhas current and warm core rings characterized by strong currents and hot temperatures where the SST is weakly correlated with SSH, (ii) the return Agulhas current with lower temperatures and currents where the SST is an active tracer, (iii) local front regions where strong SST gradients do not seem to affect the current velocities, and (iv) a weaker dynamical mode where SST is strongly correlated to SSH. Our study complements previous theoretical studies which showed that mesoscale upper ocean dynamics may be characterized by a linear coupling between SST and SSH (cf. [START_REF] Lacasce | Estimating subsurface horizontal and vertical velocities from sea-surface temperature[END_REF], [START_REF] Lapeyre | Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory[END_REF], [START_REF] Isern-Fontanet | Three-dimensional reconstruction of oceanic mesoscale currents from surface information[END_REF], [START_REF] Klein | Diagnosis of vertical velocities in the upper ocean from high resolution sea surface height[END_REF], [START_REF] Hausmann | The observed signature of mesoscale eddies in sea surface temperature and the associated heat transport[END_REF]). Following a fully observation-driven framework, the proposed latent regression model enabled us to identify different dynamical modes and to track the spacetime extension of each dynamical mode. The reported results clearly pointed out the requirement for considering a mixture model to decompose the space-time variabilities of the ocean surface dynamics. Regarding methodological aspects, it may be pointed out that EOF-based schemes (cf. [START_REF] Leuliette | Coupled pattern analysis of sea surface temperature and TOPEX/Poseidon sea surface height[END_REF], [START_REF] Casey | Sea surface temperature and sea surface height variability in the North Pacific Ocean from 1993 to 1999[END_REF]) could not reveal such non-stationary space-time variabilities. Joint EOF scheme typically decomposes a global linear mapping between the analyzed fields according to principal modes.

Regarding our future work, we will further investigate latent class regression models with additional regressors. Among others, it seems appealing to explore how time-lagged SST features and other geophysical fields such as wind speed and chlorophyll-a concentration could improve the estimation of SSH and surface currents. We also plan to apply the proposed model to other strongly active ocean regions such as the Gulf Stream system. Our objective will be to determine shared and/or system-specific dynamical modes. Future work will also investigate more detailed physical interpretation of the identified dynamical modes, especially in terms of spectral characteristics. Besides, such an improved model shall then possibly address both (i) the higher resolution prediction of mesoscale ocean surface currents from SST spatio-temporal fields (cf. [START_REF] Chen | Surface velocity estimation from satellite imagery using displaced frame central difference equation[END_REF]), and (ii) the extraction of new local and global descriptors of ocean surface dynamics from satellite sea surface observations (cf. [START_REF] Ba | Statistical descriptors of ocean regimes from the geometric regularity of SST observations[END_REF]).

Fig. 1 .

 1 Fig. 1. Satellite data of SST (a) and SSH with the associated geostrophic surface currents (b) the 1 st of September, 2004 within the Agulhas current.

Fig. 4 .

 4 Fig. 4. Distribution of the SST (a), SSH (b), surface current (U,V) norm (c) and direction (d). The results are given for each dynamical mode within the Agulhas current for the whole year 2004.

  where Y(s, t) encodes the local SSH variability through a 3dimensional vector formed by the SSH value and the surface current (U,V) at spatio-temporal location (s, t) and X(s, t) is the vectorized version of the local SST patch centered in s at time t (cf. Fig.2). It may be noted that we encoded local SSH variations at spatio-temporal location (s, t) through the surface currents which are computed as the spatial derivatives of the SSH field. As such, it constrains the method to account for spatial regularity. The p × 3 regression coefficient matrix H k (s, t) associated with dynamical mode k is corresponding to the local version of H k in Eq. (2). It corresponds to three vectorized versions of spatial convolution matrices. Here, p defines the size of the local SST neighborhood and is set according to the Rossby radius of the study region, i.e. the mean size of the mesoscale ocean structures like eddies. For the Agulhas current region, we set it up to 200 km, i.e. p = 81 for the spatial resolution of the considered data.

		X(s i , t i ) X(s i , t i )
	√ p	Y (s i , t i ) Y (s i , t i )	X(s j , t i ) X(s j , t i )
		Y (s j , t i ) Y (s j , t i )

Fig.

2

. Sketch of SST patches (contour plot), noted X, and the corresponding SSHs (scatter plot) and surface currents (quiver) noted Y at the central location s i and s j at time t i .

ACKNOWLEDGMENTS

We would like to thank the Archiving, Validation and Interpretation of Satellite Oceanographic (AVISO) and the Remote Sensing System (RSS) projects for respectively providing the altimeter-derived SSH and surface current, and the microwave SST data. This study was partially supported by the LabexMER, ANR RedHots and Institut Mines-Telecom.