Two-Particle Dispersion in Weakly Turbulent Thermal Convection

We present results from a numerical study of particle dispersion in the weakly nonlinear regime of Rayleigh-B\'enard convection of a fluid with Prandtl number around unity, where bi-stability between ideal straight convection rolls and weak turbulence in the form of Spiral Defect Chaos exists. While Lagrangian pair statistics has become a common tool for studying fully developed turbulent flows at high Reynolds numbers, we show that key characteristics of mass transport can also be found in convection systems that show no or weak turbulence. Specifically, for short times, we find Richardson-like $t^3$ scaling of pair dispersion. We explain our findings quantitatively with a model that captures the interplay of advection and diffusion. For long times we observe diffusion-like dispersion of particles that becomes independent of the individual particles' stochastic movements. The spreading rate is found to depend on the degree of spatio-temporal chaos.


Introduction
A defining characteristic of fluid flow is its ability to transport mass  [1,2]. One well-established description of transport is the dispersion of fluid elements from an initial separation r 0 , i.e. the study of their mean square separation á ñ r 2 with time. A substantive description of transport has been derived for fully developed, homogeneous, and isotropic turbulence that has a well developed inertial range  [2]. In many types of such turbulent systems  [2][3][4], dispersion has been shown to be ballistic ( ) á ñr t 2 2 or Richardson-like ( ) á ñr t 2 3  [5][6][7]. However, many flows in nature show spatio-temporal chaos or weak turbulence that do not display the scale invariance of fully developed hydrodynamic turbulence.
Here, we report numerical and analytical findings for Lagrangian dispersion of diffusing passive tracers in a vertically strongly confined flow system at very low Reynolds numbers. For the two-particle dispersion, we find scaling relationships and trace them back to simple arguments and mechanisms. The system we studied is large aspect ratio thermal convection of an incompressible fluid of Prandtl number unity, where we can clearly differentiate between effects at length scales smaller than the vertical length scale, and the influence of the vertical confinement that affects horizontal transport for larger distances.  Here, α is the volume expansion coefficient, g is the magnitude of gravitational acceleration, d the vertical distance between the plates and DT is the temperature difference between the plates (where a positive DT corresponds to the lower plate being warmer, i.e. an unstable fluid layering). Onset of convection occurs at the critical Rayleigh number » Ra 1708 c . Closely above onset, the flow forms convection rolls  [8] that, depending on initial conditions, align in parallel or display a time-dependent state known as spiral defect chaos  [9].
The Bénard system allows us to realize, in the same geometry, two qualitatively different flow patterns to study transport (figure 1). The straight roll pattern has a discrete and a continuous translational symmetry in the two horizontal directions and is time independent, which makes the flow pattern easy to describe. However, many systems in nature do not have that high a degree of regularity. The periodicity and translation invariance are broken in the state of spiral defect chaos.
For our study, we numerically solved the incompressible Navier-Stokes equations and heat equation under the Boussinesq approximation  [10] in a 2D-periodic domain of edge length = L d 100 or = L d 200 (with d the vertical system size). Straight rolls or spiral defects were selected by imposing infinitesimal temperature fluctuations to the initial conditions. The equations were discretized using a pseudospectral Galerkin scheme originally developed by Werner Pesch, and each flow component expressed at a resolution betweeń1 024 1024 2 and´2048 2048 5 spectral basis functions. Time integration used an implicit-explicit scheme of time step =t d 10 2 . Details of the GPU-based parallel numerical solver can be found in the appendix.

Dispersion of diffusing passive tracers
In this flow, we followed the trajectories of passive massless particles which, in addition to being advected with the fluid velocity u, experienced an individual Brownian diffusion with diffusion coefficient . For a small, but finite time interval t d , a particle's position x is described by the Langevin-like advection-diffusion relation where ( ) h t is a random jump direction drawn from a Gaussian distribution of standard deviation  s = t 6 d with no temporal correlation.
Particles were advected under a 4th-order Runge-Kutta time stepping scheme, where the particle velocity u i was retrieved via trilinear interpolation from sampling the spectral flow field at a resolution of1 024 1024 128, making use of the fact that the flow field ( ) u x t , changed on a time scale much larger than the advection time scale. The random diffusion jump ( ) h t was added to the particle position after each time step t d . The total number of advected particles was between 10 3 and 10 5 for each simulation run.
We considered the pair dispersion ( ) R t 2 , defined as x r ij j i , and á¼ñ denotes an averaging over an ensemble of particle pairs i j , with ( ) = r t 0 ij 0 . Particles were seeded into the flow at midheight, i.e. z=0, at random locations x y , . Note here that the presence of diffusion allows us to set the initial separation to zero, whereas nondiffusing tracers need a small initial separation for their deterministic trajectories to ever separate. We wish to point out that as an alternative measure of Lagrangian dispersion, one could take the square separation of a diffusing particle from a reference point that is purely advected by the flow. Both definitions give similar results in practice, with the single-particle dispersion about half the pair dispersion between two diffusing particles for very short and very long time scales, and the two being equal for intermediate times.
Unless noted otherwise, we measure distances in terms of the system height d and particle diffusivities in terms of the thermal diffusivity κ. As a velocity scale, we take the average flow speed in the system where á¼ñ W denotes a spatial average over the volume Ω. This velocity scale depends on the flow configuration, i.e. both the external parameters (such as the Rayleigh number), and the flow pattern. The flow velocity and system height define the vertical advection time d/U, which we take as our time scale.

Results
We find that pair dispersion in large aspect ratio Rayleigh-Bénard convection with straight rolls and spiral defects follows three different regimes of power law behavior (figure 2). For very short times, we observe normal diffusion (R t 2 ), followed by a steep t 3 -dispersion up to a plateau at For very long times, we find asymptotic behavior towards a normal diffusion, but with a much higher effective ). (b) spiral defect chaos, at Ra=2732 with = U d t 4.9 v . Data averaged from 10 independent realizations of the convective system, with 10 3 particles in each.
we propose an analytical model also known for shear flows [12] that is in excellent quantitative agreement with the observations (section 3.3). After that, we discuss the long-term dispersion for  R d (section 3.4).

Diffusive dispersion on a short time scale
For short time scales much smaller than the vertical advection time as would be the case for the relative dispersion of two points that diffuse in  3 in the absence of advection. This is unsurprising, as the flow field is nearly homogeneous on such small length scales and dispersion is a Lagrangian measure.

Dispersion on intermediate time scales
Starting from » t d U, dispersion steepens and displays the behavior The t 3 -regime starts at the same time for all particle diffusivities and in both flow configurations, and ends at dispersions around 2 . This makes particles at very low diffusivities show the scaling over more than one order of magnitude in time, while it does not appear for higher . This behavior is the same regardless of whether the flow forms straight rolls or spiral defects. If we measure time in terms of the vertical advection time d/U, we observe that the dispersion data for   2050 Ra 2732 at fixed  collapses in the range . This apparent insensitivity to the flow pattern and also the observed time scale suggests that the underlying mechanism is dominated by the convective roll structures.

Analytical model for dispersion on short and intermediate time scales
The observed features can be explained by comparing the convection system to a simple shear flow of comparable local shear rate. Dispersion in these effectively two-dimensional flows is well-understood (see e.g. Foister  [12] for a discussion of rotational and irrotational shear flow, or Ben-Naim et al [13] for unidirectional power-law shear flows), and the dispersion can be described analytically. Consider a shear flow figure 3(b)). In this velocity field, consider a particle initially at = x 0 0 at = t 0 0 . Then, the cartesian coordinates of the particle are normal distributed and have the moments (see e.g.  [12]) The distance between two such particles with identical initial positions that diffuse independently of each other is then The pair dispersion in a velocity field with a constant gradient is therefore The analytical model assumes a system that is infinitely extended, whereas particle motion in convection is confined. To observe the effects of confinement on the model, we repeated our numerical particle tracking for Couette-type shear flow between two parallel plates, periodic along the x-axis with periodicity d 2 . The resulting dispersion is shown in figure 4, where w = U d 5 was chosen as an estimate for the shear rate at the particle seed points: dispersion in the convection system collapses with both the numerical and analytical models for shear flow. Due to the appropriately chosen periodicity, which prevents particles from obtaining a horizontal distance larger than d, the simulation of confined shear flow reproduces the plateau that we observe in the Rayleigh-Bénard convection, in addition to the scaling behavior that we predict with the analytical model.

Effective diffusivity on long time scales
The transition to long time scales can be best visualized in a system with parallel straight rolls as shown in figure 1(a). Seeding a large number of diffusing particles into the system at one location allows us to study the development of their concentration with time. We examined the spatial distribution along the direction x, perpendicular to the roll orientation. The distribution assumes a Gaussian profile after a long time (  t d U 10 3 ), when the separation between the individual particles is much larger than the size of the convection rolls ( figure 5). For shorter times, the density distribution is modulated with the spatial periodicity of the convection rolls and deviates from a Gaussian. . The velocity gradient for the shear flow model is w = U d 5 . Note that the transition in scaling lies at = t d U 3 5 . Since for large t, the pair dispersion approaches a linear increase with time ( µ R t 2 ), we define the long-term effective diffusivity in the system as Setting the denominator to t 4 (instead of t 6 ) takes into account that, for very long times, dispersion is dominated by the two horizontal components of the particle separations. For ideal, straight rolls, Shraiman [11] predicts an effective diffusivity of where  is the Péclet number   º dU k , the wavenumber of the roll pattern (here » k 1), and s(k) a proportionality constant with ( ) » s 1 1.8. The so-defined effective diffusivity *  exists also in systems of spiral defect chaos. It is a function of both the Rayleigh number Ra of the flow, and the particle diffusivity  of the tracers. The dependence on the particle diffusivity is more clear at small Rayleigh numbers, whereas higher Rayleigh numbers lead to an increase of *  for all  ( figure 6(a)). In the dispersion enhancement with Ra, two factors play a role: locally, convection is faster at a higher Rayleigh number, the dependence is roughly µ -U Ra Ra c . Additionally, the pattern of the rolls changes due to the introduction of spiral defects between » Ra 2000 and » Ra 2300. A further increase in the Rayleigh number speeds up the temporal dynamics, creates more spiral cores and stronger mean currents through the system  [8,14].
In its dependence on , the effective diffusivity approaches a constant *  = const for small , and follows a power law *   µ 1 2 for large  (figure 6(b)). At low particle diffusivities, the effective particle transport is dominated by the pattern dynamics and much faster for very active spiral defect chaos, whereas at a high particle diffusivities, we find *   1 2 regardless of the underlying flow pattern. The transition point depends strongly on the Rayleigh number.
We note that this effect cannot be attributed to the effects of a generally faster advection speed U at higher Ra, but rather the additional effect of the changes in the convection pattern: the general behavior of *  remains unchanged when we normalize diffusivities by the advection velocities (i.e. express  via the Péclet number   º dU ).

Dispersion on short and intermediate time scales
We have resolved particle dispersion in spiral defect chaos for very short times. On these time scales, we find that dispersion is accelerated by the shear in the flow field. Particle dispersion scales with t 3 , which leads to a quick mixing of the diffusing particles throughout a convection roll.
Diffusion in a shear flow reproduces all characteristics of the short-term behavior of dispersion in the convective system: for very small times, we find 3D normal diffusion with   The transition to the t 3 -scaling occurs when the two terms t 12 and ( )  w t 4 3 2 3 have the same order of magnitude, i.e. around w = t 3 . This transition point only depends on the vorticity ω and is independent of the particle diffusivity. Afterwards, the t 3 -term dominates, and we find Due to the presence of top and bottom boundaries, the steep increase in dispersion is limited. When the separation distance of two particles is on the order of d, the velocity field cannot be modeled by a shear flow anymore, and the fast dispersion mechanism breaks down.
Recently, Benveniste and Drivas [15] have investigated backwards dispersion of diffusing particles in turbulent flow. Applying their calculation to the shear flow · w = u e z x (equation (7)) at viscosity ν predicts a backwards dispersion of for particles of diffusivity  n = , with a mean dissipation  nw á ñ = 2 . Since the simple shear flow is independent of viscosity, this result for backwards dispersion can be applied to any diffusivity and matches our description for the forward dispersion at intermediate times.

Dispersion on long time scales
When pair separation reaches the length scale of the system height, the transition between the rolls becomes relevant for transport, a model for which has been suggested by Shraiman  [11]. In accordance to that model, impurities in a Rayleigh-Bénard pattern of parallel, straight convection rolls have been found to spread throughout the rolls like in a normal diffusive process, albeit with an increased effective diffusion  [16,17]. For straight rolls, we find quantitative agreement with that model. An alternative dispersion mechanism is described by Camassa and Wiggins  [18], who describe how time variations in the straight roll flow cause dispersion in the absence of molecular diffusion. Both mechanisms play a role for particles in spiral defect chaos, where rolls are not stationary and the flow field is truly three-dimensional. Despite there being two large-scale dimensions available, there is apparently no large-scale advection as in two-dimensional turbulence [19] that would cause superdiffusive dispersion.
For small particle diffusivities, the effective diffusivity goes like *   1 2 , as has been observed numerically and experimentally by Gollub and Solomon [17] in a system of stationary convection rolls. When increasing the particle diffusivity, the dependence of the effective diffusivity on the particles' diffusivities vanishes, and the dispersion appears characteristic to the flow pattern itself. Our results complement the findings of Chiam et al [20], who studied the enhancement of scalar diffusion in a similar system of spiral defect chaos in terms of the Péclet number , and found a transition in scaling behavior at  » 10 2 for Ra=3074 and higher. This transition hints towards the two-fold effect of the stochastic fluid motion on particle dispersion: one part that dominates when particle diffusivity is small, and disperses particles purely by advecting them; and a second part that results from the interplay of advection and particle diffusion, which dominates at higher particle diffusivities (corresponding to lower Péclet numbers).

Conclusion
We have been able to identify and partially explain three scaling regimes for the pair dispersion that we encounter in laminar thermal convection flow. The low Rayleigh number creates a very uniform local vortex structure, such that the scaling regimes are easy to observe and characterize. The local transport mechanism is that of diffusive dispersion in a shear flow, whereas the long-term dispersion follows normal diffusion.
Due to the strong difference in the transport mechanisms before and after the transition from 3D motion to effectively 2D motion, this transition can be clearly observed in the dispersion as a plateau. If the large-scale flow pattern shows strong enough variations in time, we find that dispersion loses its dependence on the local details of stochastic motion, and depends only on the orientation pattern of the convection rolls, the correlation distance of which is small in spiral defect chaos (i.e. on the order of a few roll diameters). Despite the mean flows that appear in spiral defect chaos, long-term dispersion is diffusive in the horizontal directions.
It remains an open question how exactly the long-term dispersion can be described in terms of the pattern dynamics and correlation properties of the flow pattern. Since larger correlation lengths in the pattern require longer observations and larger system sizes, the simulation efforts are very high for these systems.
We assume a periodic domain in the horizontal directions and no-slip boundary conditions at the top and bottom plates. In the horizontal directions, the toroidal and poloidal velocity potentials ( f and g) and the temperature T can then be written as Fourier series. To account for the different boundary conditions in the vertical z-direction, the z-dependence of f T F , , and G is expressed as a series of sine functions. The zdependence of the poloidal velocity potential g is expressed as a series of Chandrasekhar functions. This spectral representation is used to construct the weak-form Galerkin solution of the equation system.
To solve for the time development of the system, we use an implicit-explicit pseudospectral scheme. We decompose the equations into their linear and nonlinear parts. The nonlinearity is evaluated explicitly and in real space. The linear part of the equation system, due to the orthogonality of the horizontal Fourier basis, can be expressed in terms of a block diagonal matrix. The vertical basis is not orthogonal, but the low Rayleigh number allows us to keep the number N z of vertical basis functions low (values between 2 and 5 were used). Both the block size of the linear matrix, and the number of Fourier transforms for the nonlinear part scale with N z . We add the nonlinearity in an explicit two-step Adams-Bashforth method and then use the linear matrix for an implicit Euler-step.
The simulation code is parallelized to run on a GPU. Unless stated otherwise, we use a resolution of 1024×1024 modes in the horizontal and two modes in the vertical direction to simulate a system of aspect ratio G = 200 with a time step of =t d 10 2 . Advection of the tracer particles is performed by sampling the velocity field in a regular grid and using trilinear interpolation to get the fluid velocity at a particle's position. In the horizontal directions, the number of sampling points is chosen equal to the number of modes, while vertically, 128 sampling points are chosen. Particles are advected in a 4th-order Runge-Kutta scheme. The final position after each time step is modified by adding a Gaussian-distributed random variable with a standard deviation of  s = t 2 d to each position component. Random variables are drawn from independent random number generators for the spatial directions and the different particles. Pseudorandom number generation uses the XORWOW algorithm.