Extreme-value statistics from Lagrangian convex hull analysis for homogeneous turbulent Boussinesq convection and MHD convection

We investigate the utility of the convex hull of many Lagrangian tracers to analyze transport properties of turbulent flows with different anisotropy. In direct numerical simulations of statistically homogeneous and stationary Navier-Stokes turbulence, neutral fluid Boussinesq convection, and MHD Boussinesq convection a comparison with Lagrangian pair dispersion shows that convex hull statistics capture the asymptotic dispersive behavior of a large group of passive tracer particles. Moreover, convex hull analysis provides additional information on the sub-ensemble of tracers that on average disperse most efficiently in the form of extreme value statistics and flow anisotropy via the geometric properties of the convex hulls. We use the convex hull surface geometry to examine the anisotropy that occurs in turbulent convection. Applying extreme value theory, we show that the maximal square extensions of convex hull vertices are well described by a classic extreme value distribution, the Gumbel distribution. During turbulent convection, intermittent convective plumes grow and accelerate the dispersion of Lagrangian tracers. Convex hull analysis yields information that supplements standard Lagrangian analysis of coherent turbulent structures and their influence on the global statistics of the flow.


I. INTRODUCTION
Turbulent transport governs the spreading of contaminants in the environment, mixing of chemical constituents in combustion engines or in stellar interiors, accretion in proto-stellar molecular clouds, acceleration of cosmic rays, and escape of hot particles from fusion machines. Because of its wide relevance, a fundamental characterization of the dispersive properties of turbulent flows is of practical interest to physicists and engineers. Here we examine the broadly relevant case of dispersion of Lagrangian tracer particles in statistically homogeneous but not necessarily isotropic turbulence.
The Lagrangian viewpoint is particularly suited to the investigation of transport in turbulent fluids. A Lagrangian description of turbulence is based on following the paths of passive tracer particles in a turbulent flow. Single-particle diffusion, as originally addressed by Taylor 1 , provides a basic characterization of a flow's transport properties 2 . A more complete characterization of the turbulent transport has conventionally been formed from the relative dispersion of two, three, or four particles [3][4][5][6][7][8][9][10][11] . However, in astrophysical environments where the effects of magnetic fields, rotation, or gravity are often significant, the more complex nature of statistically anisotropic or even inhomogenous nonlinear dynamics warrants additional examination. Dispersion in dynamically anisotropic systems such as vigorously convecting flows [12][13][14][15][16] where preferred directions exist and spatially coherent, persistent structures like convective a) j.l.pratt@exeter.ac.uk arXiv:1605.05983v4 [physics.flu-dyn] 7 May 2017 plumes can form, motivates the present consideration of a complementary diagnostic based on a different Lagrangian concept: the convex hull 17 of a n-particle group (n 4).
The convex hull is the smallest convex polygon that encloses a group of particles; two dimensional convex hulls are pictured in FIG. 1. Convex hull analysis of turbulent dispersion is similar in spirit to following a drop of dye as it spreads in a fluid, or following a puff of smoke as it spreads in the air, both classical fluid dynamics problems [18][19][20] . A large group of tracer particles can be marked, similarly to adding a drop of dye to a fluid flow, so that the same particles can be identified at all later times. Using the convex hull, a size for the group of tracer particles that were marked can be calculated at each time.
The convex hull yields statistical information about a class of Lagrangian particles that is not equivalent to preselected tracer particle groups like particle pairs or tetrads. These standard Lagrangian multi-particle statistics represent a fixed and unique structural relationship between specific tracer particles. The evolution of particle-pair structures, expressed e.g. as separation and orientation, is analyzed as the pair of particles is advected by the fluid. In contrast, the convex hull does not establish a unique link between the tracers that generate it, but continuously selects from a predefined group, based on which tracer particles have ventured furthest from the geometrical center of the ensemble. Unlike particle pairs or tetrads, the particles that constitute the convex hull are dynamically changing. The definition of the convex hull thus corresponds to a filtering based on the entire dynamical past of each particle in the group. The convex hull captures the extremes of the excursions of a group of particles, information relevant to the non-Gaussian aspects of the dynamics. The behavior of particles that do not exhibit the fastest dispersion is filtered by the convex hull, allowing a classification of particle dynamics with regard to their dispersion efficiency. In this work we begin to explore this link to extreme value theory, which has the potential to provide new physical insight for turbulent diffusion. The dynamical relation between the Lagrangian particle population forming the convex hull and the bulk ensemble of tracer particles enclosed by it represents another aspect of this diagnostic that could be exploited in investigations of turbulent structure formation.
In recent years, convex hull calculations have been used to study diverse topics such as the size of spreading GPSenabled drifters moving on the surfaces of lakes and rivers 21,22 , star-forming clusters 23 , forest fires 24 , proteins 25,26 , or clusters of contaminant particles 27 . Studies of the relationships between random walks, anomalous diffusion, extreme statistics and convex hulls have been motivated by animal home ranges [28][29][30][31][32][33][34] . Convex hulls have also been used to study analytical statistics of Burgers turbulence by analogy with Brownian motion [35][36][37] .
MHD turbulence 38,39 and turbulence during hydrodynamic convection [40][41][42] are areas where statistical analysis of Lagrangian particles has begun to be applied only recently. This work presents new Lagrangian results from threedimensional direct numerical simulations of turbulent MHD Boussinesq convection, and compares them with turbulent hydrodynamic Boussinesq convection and homogeneous isotropic turbulence. It is structured as follows. In Section II we describe the fluid simulations. In Sections III we present standard Lagrangian pair dispersion and discuss the results of these widely-used statistical tools for convective flows. In Section IV we describe the convex hull analysis that we perform on groups of many Lagrangian tracer particles. We perform several basic checks on our convex hull calculations. We then compare the dispersion curves obtained from convex hulls of large groups of particles with the expected scalings for particle-pair dispersion. In Section V we demonstrate how the convex hull can be used to examine anisotropy. In Section VI we apply extreme value theory, and show that the maximal square extensions of convex hull vertices are well described by a classic extreme value distribution, the Gumbel distribution. In Section VII we summarize the results of this validation study, and our extreme value statistics. We discuss the potential uses and benefits of convex hull analysis.

II. SIMULATIONS
We investigate three different types of turbulent systems: forced homogeneous isotropic Navier-Stokes turbulence (simulation NST) 43 , Boussinesq convection in a neutral fluid (simulation HC), and Boussinesq convection in an electrically conducting fluid (simulation MC) 12,44 . These simulations are not designed for close comparison, but produced for a broad exploration of the convex hull analysis. In each of these direct numerical simulations, the equations are solved using a pseudospectral method in a cubic simulation volume with a side of length 2π. The These equations include the solenoidal velocity field v, vorticity ω = ∇ × v, magnetic field B, and current j = ∇ × B.
The quantity θ denotes the temperature fluctuation about a linear mean temperature profile T 0 (z) where z is the direction of gravity. In eq. (3) this mean temperature gradient provides the convective drive of the system. In eq.
(1), the term including the temperature fluctuation θ is the buoyancy force. The vector g 0 is a unit vector in the direction of gravity. Three dimensionless parameters appear in the equations:ν,η, andκ. They derive from the kinematic viscosity ν, the magnetic diffusivity η, and thermal diffusivity κ. For simulation HC, the magnetic field B is set to zero. For simulation NST, both magnetic field terms and temperature terms are zero. A fixed time step and a trapezoidal leapfrog method 45 are used for the time-integration for simulation NST. The Boussinesq convection simulations HC and MC are integrated in time using a low-storage 3rdorder Runge Kutta scheme 46 and an adaptive time step, which allows for better time resolution of large fluctuations that occur during convection.
In this work we discuss turbulent dispersion in an incompressible fluid, where conservation of volume is a primitive concept. A volume of fluid that is convex at an initial time will occupy the same volume after a period of dynamic development but will generally change its shape and lose its convexity. Lagrangian tracer particles that are contained in the initial volume are marked so that they can be followed for the entire time of the simulation. At any later time, the volume of the convex hull of that group of marked particles is generally not conserved. This is illustrated in FIG. 1 for a group of particles, and for snapshots taken at three times. The growth of surface area and volume are natural concepts for convex hulls.
A summary of the fundamental parameters that describe each simulation is given in Table I. In this table, we define the Reynolds number to be Re = E is the kinetic energy, and the brackets indicate a time-average. We define the characteristic length scale L based on the largest-scale motions of the system in question. For statistically homogeneous turbulent convection the characteristic length scale is the instantaneous temperature gradient length scale L = T * /∇T 0 where T * is the root-mean-square of temperature fluctuations and ∇T 0 is the constant vertical mean temperature gradient 47 . For non-convective statistically homogeneous turbulent flows, the characteristic length scale is a dimensional estimate of the size of the largest eddies, is the time-averaged rate of kinetic energy dissipation. The magnetic Reynolds number is defined from the Reynolds number and the magnetic Prandtl number, i.e. Re m = Pr m Re. We measure length in units of the Kolmogorov microscale η kol = (ν 3 / v ) 1/4 and also make use of the Kolmogorov time-scale τ η = (ν/ v ) 1/2 . The Kolmogorov microscale multiplied by k max , the highest wavenumber in the simulation, is often used to test whether a simulation is adequately resolved on small spatial scales. In this work all of the simulations fulfill the standard criterion based on the Kolmogorov microscale (k max η kol > 1.5) for adequate spatial resolution 48 . The Reynolds numbers in Table I are on the order of 10 3 ; the Reynolds numbers and Kolmogorov microscales in Table I are in the same range as current studies of moderately turbulent flows 49,50 .
Formulation of boundary conditions for simulations of turbulent flows is delicate because boundaries strongly influence the structure and dynamics of the flow. For homogeneous isotropic turbulence, it is standard to employ boundary conditions that are periodic in x, y, and z. These fully periodic boundary conditions are used for simulation NST. For convection simulations the choice of fully periodic boundary conditions (also called homogeneous Rayleigh-Bénard boundary conditions) allows macroscopic elevator instabilities to form 51 . These instabilities destroy the natural pattern of the original turbulent flow field. The convection simulations discussed in this work use quasiperiodic rather than fully periodic boundary conditions. In quasi-periodic boundary conditions the only additional constraint is the explicit suppression of mean flows parallel to gravity, which are removed at each time step. Because our simulations are pseudospectral, the mean flow is straightforwardly isolated as the z component of the k = (0, 0, 0) mode in Fourier space, which corresponds to the volume-averaged velocity in the z-direction. Quasi-periodic boundary conditions combine the conceptual simplicity of statistical homogeneity with a physically natural convective driving of the turbulent flow. These boundary conditions do not enforce a large-scale structuring of the turbulent flow, such as the convection-cell pattern observed when Rayleigh-Bénard boundary conditions are used. In the quasi-periodic simulations presented in this work, we find no evidence of the macroscopic elevator instability although we follow the evolution of the flow for long times. Quasi-periodic boundary conditions allow for direct comparison with simulations that use fully periodic boundary conditions. In simulation NST the modes 2.5 < k < 3.5 are forced using Ornstein-Uhlenbeck processes with a finite timecorrelation on the order of the autocorrelation time of the velocity field (for further details of this forcing method, see Eswaran and Pope 52 ). The convection simulations HC and MC are Boussinesq systems driven solely by a constant temperature gradient in the vertical direction. The magnetic field present in simulation MC is generated self-consistently by the flow from a small random seed field through small-scale dynamo action. The system is evolved until a statistically stationary state is reached. For Boussinesq convection, a length scale that characterizes the scaledependent importance of convective driving is the Bolgiano-Obukhov length, bo = T , where T is the average rate of thermal energy dissipation. This length scale separates convectively-driven scales of the flow > bo from the range of scales where the temperature fluctuations behave as a passive scalar < bo . In Table I this length scale is averaged over the simulation time, normalized to the height of the simulation volume, and recorded as¯ bo . The table also includes the mean Alfvén ratio, r A = E v /E b , the time-average of the kinetic energy divided by the magnetic energy E b = B 2 /2. In the present numerical experiments, Navier-Stokes turbulence displays the weakest form of spatial coherence while Boussinesq magneto-convection exhibits anisotropy with regard to the direction of gravity as well as the occurrence of large-scale spatially-coherent structures. Additionally, a dynamical anisotropy arises because of the presence of magnetic fields.
The positions of Lagrangian tracer particles are initialized in a homogeneous random distribution at a time when the turbulent flow is in a statistically stationary steady state. The total number of particles in the simulation, n p , is listed in Table I. We use at least a million particles for a 512 3 grid. This is a standard spatial density of tracer particles used to describe homogeneous turbulence [53][54][55] . The Lagrangian statistics we produce have been tested and found to be well-resolved in space and time; we reproduce these statistics with half the particles. At each time step the particle velocities are interpolated from the instantaneous Eulerian velocity field using either a trilinear (for simulations HC and MC) or tricubic (for simulation NST) polynomial interpolation scheme. Particle positions are calculated by numerical integration of the equations of motion using a predictor-corrector method. For the convex hull calculations, the Lagrangian particle data is resampled at a rate of approximately τ η /10 for simulations NST and HC. The rate of sampling for simulation MC was smaller by a factor of 10, and this was not found to impact the dispersive results examined here. Each simulation is run for a sufficient time that Lagrangian particle pairs have separated, on average, at least by the length of the simulation volume. We call this time the Lagrangian crossing time, LCT, and it is listed in the table in units of the Kolmogorov time scale. Lagrangian particle pair dispersion statistics exhibit a diffusive trend near this time since the velocity fluctuations over this time and distance exhibit low correlation.

III. PAIR DISPERSION OF LAGRANGIAN TRACER PARTICLES DURING HOMOGENEOUS BOUSSINESQ CONVECTION
This section presents results for particle-pair dispersion during homogeneous Boussinesq convection for comparison with many-particle dispersion calculated from a convex hull analysis. For an introduction to the rich field of Lagrangian particle-pair dispersion, we refer the reader to the review of Salazar and Collins 56 . In addition to this review, several more recent works [57][58][59] propose new dispersion phenomenologies based on locally ballistic dynamics, an alternative to the classical idea of turbulent diffusion exhibiting scale-dependent diffusivity 19 . Here we briefly recall the basic argument for scaling regimes of pair dispersion. For times short compared with the autocorrelation time of the Lagrangian velocities, the relative velocity of the particles is approximately constant. The mean-squared separation of a pair of Lagrangian particles is therefore expected to grow quadratically with time for a short time. This is called the ballistic or Batchelor regime. The extent of the ballistic regime is known to depend on the initial separation of the particle pair, ∆ 0 , due to a finite correlation of ∆ 0 and the root-mean-square (RMS) velocity fluctuations on this scale v ∆0 . Recent theoretical 57-59 and experimental 60,61 works make use of a key time scale linked to ∆ 0 , the initial nonlinear turnover time τ 0 ≡ ∆ 0 /v ∆0 . In the inertial range of Navier-Stokes turbulence this initial turnover time can be estimated as τ 0 ∼ v 2 ∆0 /(2 v ). For times much larger than the autocorrelation time of the Lagrangian velocity, the velocities of a pair of Lagrangian particles are statistically independent. The mean-squared separation of a pair of Lagrangian particles is expected to grow linearly with time. This is typically called the diffusive regime. In between the ballistic regime and the diffusive regime is a period of time where the mean-squared separation of particle pairs can grow cubically with time. This is typically called the Richardson-Obukhov regime. The temporal separation of the ballistic and Richardson-Obukhov regimes 59 can be estimated by τ 0 . Achieving a clear Richardson-Obukhov regime in direct numerical simulations depends on the initial separation of particles as well as the size of the inertial range, and is the subject of current ongoing research for Navier-Stokes turbulence. For this reason, and due to the limited extent of the inertial scaling range that is expected for the Reynolds numbers we obtain, we make no claims of observing a Richardson-Obukhov regime in the present convection simulations. We compute the initial turnover time via a one-dimensional Eulerian kinetic energy spectrum as Although the moderate Reynolds numbers of the present simulations are far away from values where a true inertial range, devoid of influences from largest or smallest scales of the flows could be realized, for descriptive convenience we will apply this term to the interval of time-scales between the ballistic and the diffusive regime. FIG. 2 illustrates the Lagrangian particlepair dispersion for simulations HC and MC, both driven with homogeneous Boussinesq convection characterized by a large Bolgiano-Obukhov length. In this figure, thin solid lines indicate Batchelor scaling ∼ t 2 and diffusive scaling ∼ t, around the shortest and longest timescales, respectively. For both cases, HC and MC, we have ∆ 0 η, and thus τ 0 τ η . Indeed, both curves deviate from t 2 -scaling after τ 0 . For intermediate times 10τ 0 (t − t 0 ) 100τ 0 , they display a phase of fast separation which eventually levels off toward diffusive dispersion. The onset of fast pair separation in convection at approximately 10τ 0 is delayed compared to Navier-Stokes turbulence where it has been observed 59 to begin at (t − t 0 ) τ 0 . In a simulation of convection an anisotropy exists between the direction of the mean temperature gradient and the direction perpendicular. The separation of particle pairs evolves differently in these two directions; the separation of particle pairs can also evolve differently depending on whether the pair of particles are initially separated in the direction of the mean temperature gradient or perpendicular to it.
During Boussinesq convection with large Bolgiano-Obukhov length the Batchelor regime for pair separations looks similar to randomly forced hydrodynamic turbulence driven at the large scales, as shown by e.g. Sawford 62 , Yeung and Borgas 63 . During the diffusive regime, large-scale flow structures associated with large Bolgiano-Obukhov length Boussinesq convection clearly affect the pair dispersion curve. The dispersion curve does not look as smooth as the result obtained from randomly forced hydrodynamic turbulence driven at the large scales. This is not surprising because the separation of the particle pairs has reached sizes comparable to the large-scale convective plumes. We note that although our convection simulations use quasi-periodic boundary conditions, FIG. 2 is not qualitatively different from figure 2 of Schumacher 41 , which presents Lagrangian dispersion during Rayleigh-Bénard convection. For pair dispersion in simulations HC and MC, extensive averaging over different flow realizations would be necessary to achieve a perfectly smooth and universal result, free from the influence of intermittent plumes or large-scale magnetic structures.

A. Description of the convex hull calculations
We seed a number of tracer particles in the simulation volume, which produces a fixed density of tracer particles. In simulations NST, HC, and MC the number of tracer particles and their density is based on the number needed to produce well-resolved Lagrangian pair dispersion statistics. A convex hull analysis could potentially make use of a significantly higher density of tracer particles. To calculate a convex hull, we select and mark a group of Lagrangian tracer particles initially contained in a small cubic sub-volume of our simulation. The initial length scale of the group of particles hull is calculated as the side-length of an initial cubic sub-volume; in the limit where the group consists of only two particles, hull would be equivalent to ∆ 0 , the initial separation of a particle pair. For the density of tracer particles in simulations NST, HC, and MC, hull varies between 20 and 30 η kol . The dependence of convex hull statistics on the initial length scale and density of the particle group is examined in Appendix A.
Selection of each particle group based on the initial position of the tracer particles yields groups that contain nearly the same number of particles, with random variation of approximately 20% based on the homogeneous random initialization of the Lagrangian tracer particles. The average number of particles in a group, n pch , listed for simulations NST, HC, and MC in Table I, is between 24 and 214. We follow the N hull convex hulls (cf. Table I) of the marked particle groups for the span of the simulation. The required calculation of the hulls for each time-step is performed using the standard QuickHull algorithm 64,65 , implemented in the function convhulln in the package geometry publicly available for R, from the R Project for Statistical Computing 66,67 . The surface area and volume of the convex hulls are obtained based on a Delaunay triangulation of the hull vertices. We stop tracking the convex hull of a group of particles when the Lagrangian crossing time, LCT, is reached to avoid the possibility of numerical artifacts due to the periodicity of the simulation volume. The initial positions of particle groups could be chosen in regions of special interest in the flow, but in this work we restrict ourselves to a homogeneous initial distribution of the groups. For each simulation the ensemble of particle groups is initially selected to fill completely a horizontal slab. The total number of groups of particles that we analyze using convex hulls is listed as N hulls in Table I. This large number of convex hulls is more than are required for statistical convergence of average quantities, but allows us to capture some statistically rare flow features.
As any pair of particles separates in a turbulent flow, the particles move with the small-scale fluctuations of the velocity field. The distance between the two particles increases monotonically in time on average, but any specific pair of particles will produce an erratic, noisy signal. If a convex hull is defined by a very small group of particles, then most of the particles define the surface of the convex hull. These particles on the surface of the convex hull are called vertices of the convex hull. In the situation where most of the particles are vertices, the convex hull, like the particle-pair distance, shrinks or grows erratically as its component tracer particles move in the turbulent flow. The limit where groups contain only small numbers of particles is of little physical interest for convex hull analysis, because particle pairs or particle tetrahedra already provide useful dispersion information.
In simulations NST, HC, and MC we examine the relative dynamics of larger groups of particles. If a particle that is a vertex of the convex hull moves inward toward the center of the larger group of particles, it is unlikely that it will remain a vertex because of the requirement of convexity. It can become an interior particle of the convex hull. Other particles may continue to move away from the group, and the convex hull will typically continue to expand smoothly. The particles that constitute the group of vertices of the convex hull can be exchanged frequently. This is a distinctive concept for the convex hull because it provides a contrast with more common Lagrangian diagnostics such as particle pairs or particle tetrahedra. For statistics constructed from particle pairs or particle tetrahedra, the same particles define the size at each point in time.
The convex hull also intrinsically links a macroscopic length scale, the size of the convex hull, with the position of the convex hull's geometrical center. Over this length scale, the convex hull filters out tracer particles which disperse slower than its vertices, selecting the most efficiently dispersing members of the particle group.

B. Convex hull description of a group of tracer particles
A convex hull is defined by its vertices; these are the particles that dispersed the fastest in a given group of particles. Potentially this could decouple the convex hull from the enclosed particles in two ways. The number of vertices of the convex hull could become extremely small, or the majority of interior particles could detach from the convex hull vertices and clump somewhere in a subregion inside the hull. In this section we devise simple basic checks for these two scenarios.
If the particles contained in the convex hull do not spread throughout the space inside of the convex hull evenly as it grows, the convex hull will fail to characterize the full group of particles. We use the average difference between the geometric center of the convex hull, c vtx , and the virtual center of mass of the interior particles contained in the convex hull, c int , as an indicator of decoupling. This difference will not be zero, because the particles that make up the convex hull will never fill the space perfectly evenly. Since this difference will grow in time as the particles disperse, we compare it to a maximal extent of the convex hull at any point in time, defined by d = d 2 where d x is the extent of the convex hull projected on the x-direction, and d y and d z are defined similarly. FIG. 3(a) shows that the average difference between the centers normalized by the convex hull's maximal extent, δc = |c vtx − c int |/d . This normalized average difference between centers does not become larger than 40% during an initial phase (0.2 LCT t 0.4 LCT) and during the subsequent phase converges toward a quasi-constant level ranging between 15% and 20%, which is less than the standard deviation of the coordinates of the group of tracer particles for each simulation. In FIG. 3(a), time is given in terms of the Lagrangian crossing time, LCT.
The differences in the initial separation of the particle pairs in FIG. 2 (∆ 0 η kol ) and the mean initial length scale of the convex hulls ( hull 20 − 30η kol ) generate dispersion curves that reflect different ranges of temporal and spatial scales of the underlying turbulence. Because the observable dispersion regimes and their duration can change as a consequence of different ∆ 0 or hull , a direct comparison of both figures is difficult. In Navier-Stokes turbulence, the initial turnover time, τ 0 , has been shown 57,59 to signal the transition from the ballistic to the inertial range of dispersion, and thus to provide a reference scale of dispersion. We therefore use the initial turnover time τ 0 to normalize the dispersion of particles contained in a convex hull, with initial length scale hull . It is however not expected that this normalization can eliminate the physical differences between isotropic Navier-Stokes and anisotropic convective systems. Moreover, it is not clear whether the universality of τ 0 extends beyond the transition from ballistic to Richardson-like dispersion.
Close examination of FIGs. 2 and 3 shows that the initial phase, during which the average difference in the centers of convex hull and interior particles increases to a maximal value, extends into the fast separation regime of particle pair dispersion. The subsequent phase of decreasing δc corresponds to separation scales near to and in the diffusive regime. These signatures, as well as the sharp transients evident between phases of the evolution of δc in FIG. 3, indicate a potential utility of the convex hull for studies of the turbulent inertial range.
FIG. 4 reveals the distribution of the group of particles within the convex hull in the z-direction. Here the zdirection has been selected because it is the direction of the gravitational anisotropy in the convective cases; however for one-dimensional cuts in directions other than the z-direction, similar curves result. The ratio plotted in FIG. 4 is the standard deviation of the particle positions, σ p,z , divided by the extent of the hull in the z-direction. This ratio would be small if many of the tracer particles were to form a clump rather than spreading throughout the interior of the convex hull. For each of the three simulations we study, however, this quantity quickly comes to a plateau. After 0.1 to 0.2 LCT, i.e. the scales probed by the particles as they approach the inertial range, the ratio no longer decreases substantially.
In all simulations the average number of vertices of the convex hulls decreases only mildly with time; this decrease is on the order of 10% before the Lagrangian crossing time is reached. After the short initial phase up to τ 0 , the decrease in the number of convex hull vertices happens very gradually.
We conclude that on average in simulations NST, HC, and MC, the convex hulls and their interior particles do not detach from each other in a way that would render the concept of the convex hull inappropriate for characterizing a pre-selected group of many Lagrangian particles. Based on the measurements presented, a clear distinction can be made between the diffusive regime and the inertial range. The correlated inertial-range velocity fluctuations lead to changes in the relationship of convex hull vertices and interior particles. This trend is reversed as soon as the diffusive regime is reached, largely neutralizing the differences between interior particles and the convex hull vertices on inertial scales. This susceptibility of the convex hull to the different characteristic regimes of ballistic, inertial-range, and diffusive turbulent transport render this diagnostic attractive for future Lagrangian investigations of turbulence.
Apart from the ability of the convex hull to indicate different regimes of turbulent transport, the tests above also yield information about the dynamics of the turbulent velocity field. The average displacement shown in FIG. 3 quantifies anisotropic differences between the dynamics of the most efficiently dispersing convex hull vertices and the slower dispersing interior particles. On the spatial scale set by the convex hull, an anisotropic difference of the velocity fluctuations responsible for vertex and interior tracer transport is observable as a relative displacement of the centers of the group of interior particles and of those that define the convex hull. FIG. 3(b) which is a different representation of the data shown in FIG. 3(a)  HC have an initial difference of approximately 1%. Shifting the δc-curves of all three systems to a common initial level allows a qualitative comparison although this simple approach can not eliminate all dynamical differences caused by varying initial tracer separations.
The increase observed for δc is driven by the particles that are part of the surface of the convex hull, since they determine the geometric center of the convex hull. The relative motion of particles contained in the interior of the hull is driven by velocity fluctuations on scales smaller than the convex hull size. Particles at opposite locations on the surface of the convex hull will experience velocity differences on the scale of the convex hull and therefore tend to move more rapidly apart from each other than particles in the interior of the hull, which in turn determine the center of mass of the convex hull. Thus on time scales of (t − t 0 ) τ 0 a significant displacement between the geometric center and the center of mass of a group of particles can occur, evidenced in the rapid growth of δc. The relative displacement of the geometric center and the center of mass continues to grow at a slower rate for (t − t 0 ) > τ 0 . This can be attributed to a finite time correlation of the velocity fluctuations on the scale of the convex hull. In addition, as time evolves and the hull grows in size, particles in the interior of the convex hull will also experience increasing velocity fluctuations and thus some interior particles may become particles on the surface of the hull and -vice versa -particles on the surface of the convex hull can move into its interior due to engulfment by other particles. This process eventually leads to a decrease in the relative displacement of the geometric center and the center of mass as the diffusive regime is approached. A noticeble difference between the NST configuration and the convective systems HC and MC is the presence of a plateau for the NST case between 3τ 0 and 16τ 0 , while for HC and MC, δc continues to grow during this time. The different behavior may be caused by anisotropy in the convective flows HC and MC sustaining longer correlations in time for velocity fluctuations in preferential directions, which does not occur for the statistically isotropic Navier-Stokes case.
The quantity shown in FIG. 4 measures the diffusive character of the motion of the interior particles, rather than dynamical anisotropy. This measure exhibits a rapid transient around τ 0 from initial levels towards a first roughly constant plateau throughout the inertial range that finally approaches the asymptotic diffusive value around (t − t 0 ) 100τ 0 . Here, the inherently hydrodynamic simulations NST and HC display less variation throughout the inertial range than system MC which exhibits additional flow structuring due to the presence of magnetic field fluctuations. This brief interpretation allows for extensions, for example focusing on vertex dynamics or a detailed direction-specific analysis by introducing spatial projections of the hulls to narrow down the structure of the underlying anisotropic fluctuations. This will be subject of future work. C. Multi-particle dispersion using convex hull analysis Because ballistic and diffusive ranges for particle pair dispersion are typically discussed in terms of length squared, we employ analogous measures for a group of particles and convex hulls. This is intended to make comparison with dispersion curves as simple and direct as possible. We calculate a maximal ray r internal to a convex hull defined by a group of particles G: By definition, the particles i, j that contribute to the maximum in this definition are always vertices of the convex hull. If the group of particles densely filled a sphere, the convex hull would be the surface of the sphere, and the maximal ray would be the diameter of the sphere. For this reason the maximal ray is sometimes also called the diameter of a convex hull. However in this work we examine anisotropic systems where the convex hull of a group of particles is not typically close to spherical; we opt for the more accurate former term. The susceptibility of the maximal ray's orientation to deformations of the convex hull can be parameterized by the RMS value of the vertex distance from the hull's geometrical center normalized by the average distance to the center (averages taken over the hull vertices), Q = σ rvtx /r vtx . If Q ≈ 0, i.e. the convex hull is close to spherical, the maximal ray can change its direction by an arbitrary amount and much faster than the autocorrelation times of the underlying turbulent fluctuations would suggest. In this case, small fluctuations of the hull radius, which can occur due to uncorrelated small-scale fluctuations, will lead to rapid changes of orientation of the maximal ray. The maximal ray is highly susceptible to anisotropic deformations of the convex hull. In contrast, a significant anisotropic deformation of the convex hull, Q 1, acts like a threshold for directional variation of the maximal ray and stabilizes its orientation. This subsection focuses on quantities specific for the convex hull and their relation to classical Lagrangian mean-square pair-separation, (∆ − ∆ 0 ) 2 .
We average the square of the maximal ray over all groups of particles; the results for each simulation are shown in FIG. 5(a). This figure demonstrates that the maximal ray, although it is not tied to the same particle pair in each tracer group, asymptotically converges to a ballistic regime signature ∼ t 2 up to approximately τ 0 , and an asymptotic diffusive regime ∼ t at long times, for all systems considered. The data shown for MC does not attain the same temporal resolution as for systems NST and MC due to a larger time step but penetrates further into the diffusive regime.
Additional length-scale estimates can be obtained from taking appropriate powers of the normalized surface area, r S = (S/(4π)) 1/2 , and the volume, r V = (3V /(4π)) 1/3 , of the convex hulls. Averaging the length scales produced by many different particle groups reveals dispersive behaviors that also tend to obey the ballistic and diffusive scaling laws. A comparison of dispersion curves produced from the surface area and volume of simulation NST are shown in FIG. 5(b). Similar to Lagrangian pair dispersion, the expected asymptotic scaling laws for ballistic and diffusive regimes are approached by the surface-and volume-based distance approximation. However, they hold over a shorter period of time than those shown in FIG. 2. Although FIG. 5(b) shows dispersion curves only for simulation NST, similar results are found for simulations HC and MC. A Richardson-Obukhov-like regime is not observed. Because achieving a clear Richardson-Obukhov regime in direct numerical simulations depends on the initial separation of particles as well as the size of the inertial range, a Richardson-Obukhov regime is not expected in our simulations. Particle filtering, an inherent property of the selection criterion of convex hull vertices, may also contribute to the lack of a clear Richardson-Obukhov regime resulting from convex hull analysis of dispersion. During early dispersion, the vertices of a convex hull tend to be particles that move away from the center of the hull most rapidly in the direction radially outward from the center of the particle group; this may explain the quasi-ballistic signature before approximately 16τ 0 (cf. FIG. 3(b)). As noted by Bianchi et al. 50 , although there is a conceptual connection between many-particle groups and particle pairs, many-particle groups provide different information when measuring dispersion scalings.
There is a fundamental difference between the maximal ray and the surface-or volume-based length approximations that becomes particularly important with regard to deformations of the convex hull: the maximal ray by definition runs along the direction of maximum extent of the convex hull. In contrast, the other two quantities yield averaged and isotropized approximations of the length scale probed by the hull, i.e. the radius of a reference sphere of same surface or volume. Spherical geometry is a natural first-order approximation of a convex hull or, more precisely, the convex polyhedron, which we use as its numerical representation, since convexity implies that the hull has no corner pointing inwards. This constraint severely restricts the complexity of the hull's surface structure, since any such corner vertex would turn into an interior point enclosed by the hull. This results in an object which can mainly be deformed by flattening of the inscribed spheroid along some direction perpendicular to the maximal ray. The convex hull is not material and therefore is not constrained by volume conservation in incompressible flow. Although the possible length definitions do not show large qualitative differences compared to the maximal ray, their behaviour relative to each other reflects the different responses of hull area and volume to deformations of the convex hull. This will be exploited in Section V.

V. RESULTS: ANISOTROPIC DYNAMICS OF CONVEX HULL VERTICES
The relationship between the surface area, S, and volume, V , of a convex hull reveals the anisotropy of vertex transport in a turbulent flow, which is of particular interest during convection, i.e. in the presence of coherent velocity structures. We introduce the non-dimensional ratio S/V 2/3 as a direct way to quantify anisotropy. Because a sphere minimizes the amount of surface area for a given volume, an absolute lower bound of 4π/(4/3π) 2/3 ≈ 4.8 exists for this non-dimensional surface-volume ratio. An anisotropic convex hull, e.g. a cigar-shaped or a pancake-shaped hull, will have a higher surface to volume ratio, so the ratio gives an impression of how non-spherical the current state of the hull is. The ratio cannot differentiate between prolate (cigar-shaped) and oblate (pancake-shaped) convex hulls, because it approaches infinity in the limit both of zero pancake thickness and infinite cigar length. Higher values indicate a basic level of anisotropic deformation. FIG. 6(a) shows the time evolution of the surface-volume ratio, averaged over all convex hulls in each simulation. Because the particle groups consist of small numbers of particles which are randomly distributed, they are not initially perfectly isotropic and do not evenly fill the cubic initial volumes; the resulting convex hulls do not form either perfect cubes or perfect spheres. Thus the surface-volume ratio initially exhibits an average value of approximately 5.6, a low value that lies between the values for perfectly spherical and perfectly cubical volumes. In all simulations, the surface-volume ratio begins to increase around t = τ η , indicating that the convex hulls typically become stretched, i.e. anisotropic, as their particles start to disperse due to turbulent fluctuations. In the Navier-Stokes case (NST) no global anisotropy exists in the flow. As expected, the average surface-volume ratio remains relatively low throughout the simulation reaching its maximal value around 10 τ η . At long times, the average surface-volume ratio returns to approximately its initial value as uncorrelated particle motion begins to eliminate anisotropic deformations of the convex hull. The changes in the surface-volume ratio also slow and it approaches a flat regime related to the diffusive trend observed in FIG. 5 at long times.
In the case of hydrodynamic Boussinesq convection (HC), the mean temperature gradient introduces a preferential direction. We would thus straightforwardly expect higher stretching of the hulls in this direction. However, FIG. 6 shows that this does not take place for the convex hulls we followed; for times greater than τ η only a slight increase occurs followed by a plateau phase up to 30τ η . Subsequently, the average surface-volume ratio quickly decreases below its initial value. The scale of the convective plumes in simulation HC are large and diffuse, as reflected by the large Bolgiano-Obukhov length bo , the smallest scale on which the cascade of thermal fluctuations is driven by buoyancy 68 . This large Bolgiano-Obukhov length indicates that smaller-scale turbulent dynamics are not driven by the anisotropic influence of buoyancy. The convex hulls do not tend to become strongly anisotropic, because the length scale of the anisotropic convection differs considerably from the scale of the convex hulls examined. The Reynolds number of HC which is less by approximately 50% as compared to the value of the NST system also explains why the HC simulation exhibits the lowest level of convex-hull anisotropy.
A different behavior is observed for the magnetohydrodynamic convection (MC) simulation since larger-scale mag- netic fluctuations have a strong impact on small-scale dynamics (in contrast to a large-scale velocity there exists no frame of reference which eliminates the magnetic field); consequently, far higher surface-volume ratios are attained than in the other two cases. In this simulation the large-scale magnetic field fluctuations result in strong local anisotropy of the small-scale velocity fluctuations [69][70][71][72][73][74][75][76] ; the consequence is considerable stretching of the convex hulls. The mean alone does not characterize the full information that the convex hull analysis can provide about anisotropy in each simulation. The shape of the probability distribution of the surface-volume ratio yields a more comprehensive picture. If all convex hulls in a simulation were perfect spheres, the distribution of the surface-volume ratio would be a delta function. However, the distributions show a strong dependence on the type of turbulence indicated by the values of distribution mean, µ, and standard deviation, σ, given in the caption of FIG. 6. In the hydrodynamic convection case the distribution is the narrowest with the lowest mean indicating the highest level of anisotropy, followed by NST and MC. The significant hull anisotropy observed for system MC is a clear indication of the additional anisotropy imposed by the slowly evolving large-scale magnetic field fluctuations on the smaller-scale velocity fluctuations. These results are not surprising and consistent with the data given in FIG. 6(a). In addition, FIG. 6(b) shows the centered and normalized distributions of the surface-volume ratio for each of our three simulations after each set of convex hulls has evolved for 20 τ η . All distributions collapse on a positively skewed functional shape suggesting a general characteristic of convex hull deformation common to all three turbulent systems.
The surface-volume ratio varies spatially in each simulation. The time evolution of this ratio for a single convex hull in simulation MC is illustrated in FIG. 7. At early times, the surface-volume ratio for this individual hull grows to considerably exceed the mean, indicating that this hull is more stretched than the average convex hull of this ensemble. This surface-volume ratio also exhibits rapid changes in time. For example, during the period between approximately 5 τ η and 10 τ η this hull goes from a more anisotropic form than average to a considerably less anisotropic form.
In FIG. 7, the surface-volume ratio is also shown as a contour plot for the set of convex hulls that fill a horizontal slab of simulation MC. Dark areas represent regions where convex hulls have grown with significant anisotropy. High spatial intermittency is also noticeable, with areas of large anisotropy bordering areas that grow more isotropically. This pattern of anisotropy remains similar over a long period of time, reflecting the strong influence of the initial configuration of the flow on local dispersion. Although we examine a small number of simulations, the non-dimensional surface-volume ratio that we introduce is clearly capable of revealing aspects of local anisotropy in turbulent flows.

VI. RESULTS: EXTREME-VALUE STATISTICS OF TURBULENT PARTICLE DISPERSION
The vertices of a convex hull are the particles that disperse fastest among a given group of particles, and the maximal ray defines a maximal dispersion of all particle pairs within the group. Thus the use of the convex hull evokes concepts from extreme value theory 77,78 . The most widely encountered distribution in extreme value theory, the Gumbel distribution 79,80 , has been frequently employed for climate modeling, including extreme rainfall and flooding [81][82][83][84][85] , extreme winds 86 , avalanches 87 , and earthquakes 88 . The Gumbel distribution has also been found to reasonably characterize the density fluctuations within galaxies [89][90][91] and in certain areas of tokamaks [92][93][94] , binding energies in liquids 95 , as well as turbulent fluctuations 96,97 . The cumulative distribution function F for the Gumbel case has the well-known form: where the location parameter µ gives the mode of the distribution, β is commonly called the scale parameter, and the median of the distribution is µ − β ln(ln (2)). Because extreme value theory typically develops as an asymptotic theory for sample sizes n ∼ ∞, convex hulls with large numbers n of particles facilitate the exploitation of extreme value theory results. We examine the square-length of the maximal ray with extreme value theory, and this choice is crucial. The square-length of the maximal ray is a fundamental scalar commonly associated with dispersion, and thus the most natural physical quantity to consider. The square-length of the maximal ray is also consistent with a simple model of Gaussian displacements. No rigid upper limit exists for the square-length of the maximal ray, and thus the Gumbel distribution is the case that would be anticipated from extreme value theory.
Because Lagrangian tracer particles move in a flow with a finite correlation in space and time, their motions are not independent. The number of particles in each group is also limited in these numerical experiments. Despite these limitations, we find that the shape of the cumulative distribution function of the square of the maximal ray is suggestive of a Gumbel distribution. This observation holds at each point in time, regardless of whether the particle groups sampled are in the ballistic regime, diffusive regime, or a transitional period of dispersion. A Gumbel distribution describes the results well, regardless of the initial length scale of the convex hull, and the initial density of particles, for the range 4η kol < hull < 64η kol that we have tested (see Appendix A). This suggests that the Gumbel distribution might provide an effective description of the probability of extremes of turbulent dispersion. The location and scale parameters can be different for different hull , and at different times in the dispersion process, although a Gumbel distribution is recovered at each time.
In addition, we consider a cumulative distribution function constructed from data at all times throughout the evolution of the convex hulls, as shown in FIG. 8. Using data from all times is a reasonable choice that produces a single form of the cumulative distribution function relevant to the entire simulation. From the perspective of the simple model of Gaussian displacement, noted above, that pragmatic choice actualizes a distribution of values of the scale parameter. Such a possibility is well known in related, but physically distinct, studies of turbulence 98 . FIG. 8(a) shows that the distribution of square-length of the maximal ray is fit well with a Gumbel distribution when physically distinct directions, perpendicular and parallel to gravity, are considered individually in the magnetoconvection simulation MC. r 2 /η kol 2 ( 10 7 ) ln[−ln[F(r 2 /η kol 2 ))]] q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q (a) ))]] q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q We found in Section V that the convex hulls in simulation MC become highly anisotropic on average. Thus the fact that a Gumbel distribution with different location and scale parameters accurately describes the extremes of dispersion in both of these physically distinct directions is a new and significant physical observation. In FIG. 8(a), we observe an ordering between the scale parameter obtained for the direction perpendicular to gravity and the direction parallel to gravity; the value of the scale parameter is larger in the direction parallel to gravity. FIG. 8(b) compares the cumulative distribution functions of the square-length of the full maximal ray of the convex hull in each simulation, and again they demonstrate the linear behavior expected of ln (− ln(F )) for the Gumbel distribution. When the ln (− ln(F )) is fit using linear regression, the value of the scale parameters are: β HC = 0.17, β NST = 0.40, and β MC = 0.65. In Section V, ordering these simulations according to the least anisotropic to the most anisotropic simulation produced: HC, NST, MC. We thus conjecture from the results in FIG. 8(a) and (b) that faster dispersion linked to anisotropy will lead to a higher value of the scale parameter.

VII. DISCUSSION
We have shown that the convex hull can be used to characterize many-particle dispersion in turbulent flows, and can reproduce scalings similar to particle-pair, and other multi-particle Lagrangian statistics. The convex hull allows us to extract dispersion behaviors that produce clear scalings from groups of tracer particles that are significantly larger than have been typically examined by multi-particle statistics. We have examined particle dispersion using convex hulls across three types of physically distinct turbulence simulations, including Navier-Stokes turbulence, Boussinesq convection, and MHD Boussinesq convection. In each of the simulations that we consider, we have shown that the convex hull describes well the dynamics of the entire group of particles. In addition, these tests yield further information about the turbulent velocity field by quantifying the dynamical differences between interior particles and convex hull vertices. Dispersion curves produced using the maximal ray of the convex hull, the surface area of the convex hull, and the volume of the convex hull produce ballistic and diffusive scalings, which can be compared with particle-pair dispersion curves. Although the convex hull has been used to calculate volumes occupied by particles in some specialized contexts 21,27 , this is the first time that the convex hull of the positions of Lagrangian tracer particles has been used as a fundamental diagnostic to obtain Lagrangian statistics of multi-particle dispersion in homogeneous turbulent flows.
In addition, we have explored the convex hull's fundamental link to extreme value statistics. We have discussed that the convex hull provides new information about extremes of dispersion that standard multi-particle statistics cannot. Convex hulls calculated from large numbers of particles provide an ideal application for extreme value theory, an asymptotic theory for large samples. Predictions based on extreme value theory are of practical use for studies of contaminants or of energetic particles, where questions about maximal dispersion are critical. Experimentally it may be simpler to track the convex hull of a large number of particles than to track all the particles in the group individually. We show that the distribution of the square length of the maximal ray of the convex hull is the Gumbel case of generalized extreme value distributions. In addition we show that for a system that is anisotropic because of MHD convection, the maximal ray in each physically distinct direction is described well by the Gumbel distribution. Because the Gumbel distribution has been successful in predicting avalanches, extreme rainfall, and extreme winds, this nontrivial new observation will provide new physical intuition for modeling anomalous dispersion.
In a second application of the convex hull analysis, we exploit the relationship between convex hull surface area and volume to examine the degree of anisotropy present in a turbulent convective flow. Our results reveal the extent of spatial variation of anisotropy. Moreover, this quantity also exhibits a probability distribution that has the same universal shape for all three considered physical systems. Convex hull analysis can easily isolate dispersive characteristics in any local region of interest, for example a region where a magnetic structure, or strong convective plume is present. Used in this way, they provide a versatile supplement to standard Lagrangian multi-particle statistics in complex turbulent flows. Because of these advantages, further investigation of the convex hull to analyze manyparticle turbulent dispersion is justified. FIG. 9. Evolution of the difference between the geometric center of the convex hull and the center of mass of the group of tracer particles, normalized by the initial length scale of the convex hull |cvtx − c int |/ hull for (a) groups of particles with four initial sizes and particle density ρ hi , and (b) groups of particles with eight initial sizes and particle density ρ low . The initial size hull is labeled in units of η kol . A thin solid line with slope 1 is indicated during time scales associated with the inertial range.
Aside from these diagnostics, we consider the growth of the difference between the geometric center of the convex hull and the center of mass of the group of tracer particles, normalized by the initial length scale of the convex hull |c vtx − c int |/ hull . This is not useful for examining whether the convex hull describes the group of particles well, because unlike the maximal extent d, hull describes the initial state of the particle group and does not change in time; however the difference in these centers provides a new quantity linked to the dispersion of many tracer particles. The evolution of this quantity is shown in FIG. 9. In the figure, the magnitude of the difference in the centers is clearly linked to hull , with larger hull leading to smaller values of |c vtx − c int |/ hull throughout dispersion. However the shape of the evolution curves for the difference in the centers appears to be approximately independent of hull . During the inertial range of time scales, the difference in centers grows linearly with time. Explaining this interesting scaling result will be the subject of future work.  10. Evolution of the mean-square maximal ray r for groups of particles with (a) four initial sizes and particle density ρ hi , and (b) eight initial sizes and particle density ρ low . The initial size hull is labeled in units of η kol . Thin solid lines with slope 2 and slope 1 indicate the ballistic and diffusive regimes respectively.
The evolution of the maximal ray r introduced in Section IV C universally exhibits a clear diffusive regime for all hull and densities tested. This is illustrated by FIG. 10. We do not expect perfect ballistic scaling of the maximal ray, because unlike a particle pair, the particles that determine the maximal ray of a convex hull can be exchanged. Despite this, we do observe a scaling reminiscent of ballistic behavior for all hull and densities tested; this likely indicates that vertex exchange is not a dominant effect during this early regime of dispersion. The slope of the dispersion curves during transitional regimes between ballistic and diffusive appears to be dependent on the initial length scale hull . This is unsurprising because it is a well-known result for particle-pairs.