Two-dimensional turbulence dispersion in a closed domain: Influence of confinement and geometry

The spreading of passive particles immersed in a two-dimensional turbulent flow confined within a closed domain is studied analytically and numerically. The primary goal is to investigate the effect of the confinement and the geometry of the container on one and two-particle Lagrangian statistics (absolute dispersion from point sources and relative dispersion of pairs of particles, respectively). The influence of the flow confinement is analysed by performing numerical experiments with numerous particles in square boxes with different sizes. The results examine the modification of the time-dependent, dispersion curves as the particles spread out (in comparison to the turbulent regimes for unbounded flows). At long times, such curves asymptotically reach a constant value of saturation as the particles fill the container. Theoretical saturation values are calculated, and the obtained formulae are tested with the numerical results. To study the influence of the domain shape, saturation values are computed analytically for different geometries (rectangles, triangles, and ellipses). To our knowledge, the obtained expressions are new. The saturation values depend on the characteristic lengths of the domain for both regular and irregular shapes. Ranges of saturated values for the different geometries are provided. The results are compared with well-known asymptotic values for unbounded flows, thus determining the influence of the closed boundaries on particle dispersion. © 2019 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5081848 The transport and fate of material or fluid properties are determined by a wide range of motion scales encompassed by turbulence. Lagrangian statistics, based on the average behaviour of a large number of discrete tracers, are a useful tool for describing turbulent transport.1 This method consists of releasing a large number of particles in the flow and recording their positions and velocities; then, the information of individual trajectories is used to determine statistical quantities that represent the overall dispersive properties of the turbulent flow. Two pioneering works on this topic date back to the 1920s. First, Taylor2 examined statistically the average separation of particles from their original position (the so-called absolute dispersion). Taylor’s study was one of the first successful attempts to develop a statistical framework for turbulent dispersion. The theory was able to explain both, the early ballistic regime (quadratic in time), and the random-walk behaviour when particle motions are decorrelated. Later, Richardson3 considered a large number of particle pairs and studied the average relative separations (relative dispersion). This approach provides information on the evolution of a cloud of particles, regardless of their initial positions. Richardson’s aim was to explain the different diffusivity (rate of dispersion) observed in flows at different length scales. In homogeneous two-dimensional (2D) turbulence, the energy spectrum has two inertial ranges separated by the forcing scale.4 A range with upscale energy flux towards larger scales (inverse energy cascade), and another one with downscale enstrophy flux down to dissipation scales (direct enstrophy cascade). Numerous authors have examined the problem of particle dispersion under the 2D dynamics. The studies of Babiano et al.5,6 (and references therein) are focused in the dispersion of individual particles and pairs, and their relationship with the energy spectrum. Provenzale7 examined the turbulent transport by coherent vortices. Experimental works with electromagnetically-forced 2D turbulence analyzed pair dispersion in the inverse energy8 and direct enstrophy9 cascades. The first case was also studied by means of AIP Advances 9, 035302 (2019); doi: 10.1063/1.5081848 9, 035302-1

The transport and fate of material or fluid properties are determined by a wide range of motion scales encompassed by turbulence. Lagrangian statistics, based on the average behaviour of a large number of discrete tracers, are a useful tool for describing turbulent transport. 1 This method consists of releasing a large number of particles in the flow and recording their positions and velocities; then, the information of individual trajectories is used to determine statistical quantities that represent the overall dispersive properties of the turbulent flow. Two pioneering works on this topic date back to the 1920s. First, Taylor 2 examined statistically the average separation of particles from their original position (the so-called absolute dispersion). Taylor's study was one of the first successful attempts to develop a statistical framework for turbulent dispersion. The theory was able to explain both, the early ballistic regime (quadratic in time), and the random-walk behaviour when particle motions are decorrelated. Later, Richardson 3 considered a large number of particle pairs and studied the average relative separations (relative dispersion). This approach provides information on the evolution of a cloud of particles, regardless of their initial positions. Richardson's aim was to explain the different diffusivity (rate of dispersion) observed in flows at different length scales.
In homogeneous two-dimensional (2D) turbulence, the energy spectrum has two inertial ranges separated by the forcing scale. 4 A range with upscale energy flux towards larger scales (inverse energy cascade), and another one with downscale enstrophy flux down to dissipation scales (direct enstrophy cascade). Numerous authors have examined the problem of particle dispersion under the 2D dynamics. The studies of Babiano et al. 5,6 (and references therein) are focused in the dispersion of individual particles and pairs, and their relationship with the energy spectrum. Provenzale 7 examined the turbulent transport by coherent vortices. Experimental works with electromagnetically-forced 2D turbulence analyzed pair dispersion in the inverse energy 8 and direct enstrophy 9 cascades. The first case was also studied by means of numerical simulations. 10 Recently, Salazar and Collins 11 devoted their review to laboratory experiments and numerical simulations on two-particle dispersion (in the context of 2D and 3D turbulent flows). From a geophysical point of view, Lagrangian statistical measures of oceanic observations are summarized in the work of LaCasce. 12 Particle dispersion, and therefore the Lagrangian statistics, can be modified by the effects of boundaries. For instance, solid walls may play an important role in the evolution of 2D turbulence since they act as sources of vorticity that is advected to the flow interior in the form of filaments and vortices. 13 As a consequence, the particle statistics are altered because the characteristic scales of motion and the size of the domain are not well separated. 14 In geophysical situations, the presence of lateral boundaries may affect the dispersive patterns of surface drifters in ocean basins, as has been observed in the Mediterranean Sea 15 and the Gulf of Mexico. 16 Another factor is the shape of the flow domain, because introduces anisotropic effects on particle dispersion, as shown in elongated basins such as the Adriatic Sea 17 and the Gulf of California. 18 In this paper, we present the Lagrangian statistics of a forced, 2D turbulent flow bounded by solid walls measured from numerical simulations and analytical considerations. The study is focused on the effects of both, the confinement (size) and geometry of the domain on time-dependent particle statistics. After describing the methods to simulate a turbulent flow and the design of dispersion experiments in section I, we present the results following two main approaches.
First, we numerically calculate one and two-particle statistics in different square boxes. In section II, absolute dispersion from an arbitrary source is shown, whilst the relative dispersion and kurtosis of pair separations are presented in section III. We will show how these statistical measures are modified because the particle trajectories are limited by the presence of the boundaries. Furthermore, as the particles fill the domain, the dispersion curves asymptotically approach a fixed value (namely, the saturated absolute dispersion, relative dispersion and kurtosis of saturation). We derive analytical formulas for the saturation values, which depend on the domain size.
Second, we describe the effect of the domain geometry in section IV. In order to do so, the saturation values are calculated for simple regular and irregular shapes. From the derived formulas we determine the ranges of saturated values for the different geometries, and compare with asymptotic values for unbounded flows. Finally, the main conclusions are presented in section V.

A. Turbulent flow simulations
We consider a 2D flow in a Cartesian coordinate system (x, y). The flow dynamics is described by the vorticity-stream function formulation of the 2D Navier-Stokes equation: where t is the time, ν the kinematic viscosity, λ is a linear friction coefficient, ω = ∂v/∂x − ∂u/∂y is the vorticity, with (u, v) the velocity components in the (x, y) directions, respectively, and F is the curl of an external forcing. The flow is assumed incompressible, so that the stream function ψ is defined in terms of the velocity field as u = ∂ψ/∂y, v = −∂ψ/∂x. The operators ∇ 2 = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 and J(a, b) = (∂a/∂x)(∂b/∂y) − (∂b/∂y)(∂a/∂x) are the Laplacian and Jacobian operators, respectively. The vorticity and stream function are related through the Poisson equation ω = −∇ 2 ψ. The flow domain is a square with side length L and no-slip boundary conditions. The fluid is set in motion by the external forcing. In this work, we consider the function F with magnitude F 0 and a chess board shape: where m, n are positive integers. The product of sine functions initially generates an array of 2m × 2n vortices (in the x and y-direction, respectively) with alternating vorticity. Moreover, it does not impose any net torque in the fluid inside the square domain. 19 The forcing term is multiplied by two auxiliary functions: h(t) initially sets the forcing equal to zero and gradually increases the magnitude of the forcing at a time scale τF. The function g(x, y) attenuates the forcing in the vicinity of the walls, thus avoiding the advection (and the consequent retention) of particles against the no-slip boundaries. The external forcing continuously injects energy into the flow. Most of the energy being transferred to larger scales (due to the inverse cascade) is dissipated by the linear damping term, thus avoiding the pile-up of energy at large scale motions. Linear friction is often used to model 3D effects, such as the influence of an underlying fluid layer, or Ekman bottom friction in rotating systems. 20 The vorticity equation (1) was solved by using a finite difference model initially developed to study the evolution of 2D flows, 21 and later on used to simulate quasi-2D turbulence with variable topography 22,23 and forced shear flows. 24 The main features of the numerical code are: the viscous terms are discretized by a centered, second-order scheme and the advective terms by an Arakawa method. The non-linear and viscous terms are discretized in time by an explicit third-order Runge-Kutta method and an implicit scheme, respectively. 21,25 To examine the influence of the domain size on particle dispersion, we performed three long-term, canonical runs with different L values and essentially the same turbulent flow, as described below. The domains are denominated as large (L = 1 in arbitrary units), medium (L = 0.5) and small (L = 0.25).
To describe the turbulent flow, consider first the large domain L = 1. The square is discretized with N = 513 2 grid points, i.e. a spatial resolution of ∆x = ∆y ≃ 0.0019. The frictional parameters are ν = 10 −6 and λ = 1.7 × 10 −3 , and the forcing ones are F 0 = 0.17, m, n = 20, δF = 0.45 and τF = 10. Note that integers m and n determine the forcing length scale, which is of order L/m ∼ 0.05. The total time of the simulation is 5000 time units with a time step of ∆t = 0.1. In all cases, we used no-slip boundary conditions. The fluid is initially at rest and then set in motion by the gradually-increasing forcing. The main parameters ν, λ, and F 0 , were chosen to generate a statistically steady turbulent flow once the forcing reaches a maximum value. Figure 1  evolution of the global energy, E = 1/2∬(u 2 + v 2 )dxdy, and enstrophy, Z = 1/2∬ω 2 dxdy. At early times, t ∼ O(τF), these quantities grow up as the forcing gradually acts. Afterwards, both of them decay and evolve around a mean value, Em = 1.496 × 10 −5 and Zm = 0.128, respectively. The corresponding fluctuations over the interval 500 ≤ t ≤ 5000 are small (about 3 and 6%, respectively). During this period the flow is considered to be in a statistically stationary state, adequate to perform experiments on particle dispersion. A typical vorticity field of the statistically stationary turbulence is shown in Figure 2a. Figure 2b shows the energy spectrum Ê(k), where k is the wave number, calculated at different times indicated with arrows in Figure 1. The spectra present a similar shape at all times, which confirms the stationarity of the turbulence. More importantly, two inertial ranges are observed, as predicted by homogeneous, 2D turbulence theory for unbounded flows. 4 These ranges are separated at the cut-off wave number k f = 131.946 (vertical solid line) determined by the external forcing. For 40 < k < k f (inertial large scales) the spectrum shows a good agreement with the theoretically predicted slope of − 5/3. The linear friction dissipates the energy that accumulates at scales of the order of the domain size (k < 40). For k f < k < 500 (inertial short scales) the spectrum shows a − 6 slope, which is steeper than the theoretical value − 3. A possible reason for this discrepancy might be the confinement of the flow in a closed domain, as reported in laboratory experiments of electromagnetically-forced quasi-2D turbulence, 26 where the flow motion is unavoidably bounded. On the other hand, some studies 27,28 pointed out that the inclusion of the linear friction term affects the enstrophy cascade causing steeper power laws (<− 3) in the spectrum that differ from the Kraichnan prediction.
A macroscopic Reynolds number defined in terms of the domain size is 29 of order urmsL/ν ≈ 3500, with the root-meansquare velocity urms = 3.57 × 10 −3 . To compare the turbulent flow developed in the simulations with medium (L = 0.5) and small (L = 0.25) domains, it is convenient to define the Reynolds number based on the forcing scale r f = 2π/k f = 0.0476. Some of the flow parameters used for the three domains are presented in Table I. Since r f ∼ L/m is the same in all cases, then the turbulent flows have a similar Re. Note that the spatial resolution is kept constant. The turbulent layer near the boundaries, which is of the order of the forcing scale, contains approximately r f /∆x ∼ 26 grid points. A classical viscous boundary layer is about ∼ 5r f /Re 1/2 , which is about 0.0192. Thus, there are at least ten grid points to solve these scales.

B. Dispersion experiments
Given the domain size (large, medium or small), we performed 20 simulations whose initial conditions are the vorticity fields in  the corresponding canonical simulation calculated at times indicated in Figure 1a (i.e. the times at which the energy spectra were calculated). The time between each vorticity field is 200, which is much larger than the Lagrangian time scale ≈ 9 (calculated in preliminary particle experiments). Using this procedure we obtained 20 independent realizations of the flow to calculate ensemble averages. For each experiment we released a set of N = 1024 particles on a square array of M × M particles (with M = 32) separated each other a distance ∆ = 0.0021 ( Figure 3). In most of the experiments, the array is placed at the centre of the domain, except in some cases where it was located near one of the boundaries. The design of the initial array of particles considers two relevant situations. First, the minimum separation between particles is smaller than the forcing scale ∆ < r f , which implies that the dispersive effects of small scales in the direct-enstrophy range are taken into account. This is particularly relevant for pair dispersion. Second, the initial array of particles is somewhat larger than the forcing scale, 31∆ > r f , which ensures that some particles are advected by different size eddies at the beginning of the simulations (avoiding the initial advection of the whole set of particles as a single blob). The position of a particle k located at (x k , y k ) is obtained by solving with k = 1, 2, . . ., N. We integrate equations (3) for a set of N particles using a second-order Runge-Kutta method with the same time step of the flow simulations. The positions of the particles are recorded every 5∆t. Therefore the statistical measures are obtained by using an ensemble of 20N = 20480 particles with a spatial resolution of t = 0.5. All the procedures described above were repeated for the small and medium-size square containers.

A. Saturated absolute dispersion
Absolute dispersion for a group of N particles is defined as the mean square separation from their initial position:

ARTICLE
scitation.org/journal/adv where ⟨⋅⟩ represents an ensemble average and t 0 is the time at which the particles are released. 7 In homogeneous turbulence, absolute dispersion is a growing function of time. For unbounded domains, this quantity behaves differently at small and large asymptotic times: where EL is the mean Lagrangian kinetic energy and D is a constant diffusion coefficient. 2,7,30 The quadratic dispersion at small times is known as ballistic regime, and the linear dispersion for longer times corresponds to the standard diffusion (or random-walk) regime. When the domain is finite, the growth of the absolute dispersion will reach a limit as the particles fill the basin. Figure 4 illustrates the particle positions calculated at two different times in a single simulation with the large domain (L = 1). At time t = 20, shortly after being released, the particles spread out from the centre of the box, and at t = 500 they saturate the domain. In this section, we calculate the limit value of absolute dispersion and test it with numerical simulations.
The saturated absolute dispersion is calculated by averaging the squared distance to the origin, a 2 = x 2 + y 2 , of an arbitrarily large number of particles filling the container. This average is obtained by integrating over the whole square domain: Thus, the saturated absolute dispersion is 1/6 of the area of an arbitrary square. This procedure is generalized in section IV for different geometries. Figure 5 shows the absolute dispersion calculated from definition (4) for the three square domains. The theoretical saturation values are indicated with horizontal dashed lines. At short times (t < 10), all curves behave similarly to the ballistic regime. During this period, the spreading of particles only depends on the turbulent flow and is not affected by the walls. At longer times, the dispersion curves slow-down, follow the standard dispersion regime for a short period, and then become uniform in time. The transition is longer for the large square because the particles take more time to saturate the domain. For the small-size square, in contrast, the transition from the ballistic regime to saturated dispersion is very fast because of its limited area. In all cases, the dispersion reaches the expected value, L 2 /6, predicted in (6).

B. Absolute dispersion from an arbitrary source
The saturated absolute dispersion was calculated for particles initially concentrated at the center of the domain (in the simulations, the initial array of particles is sufficiently small to consider that they are released from a point source at the origin). The saturation value can be calculated when the particles are released from an arbitrary location in the domain. If the source of particles is located at (x 0 , y 0 ), then the saturated absolute dispersion is Thus, the saturation value also depends on the initial location of the particles. The minimum saturation value is obtained when the particles are released from the origin, whilst the maximum value corresponds to the initial position located at one of the corners of the square (x 0 , y 0 ) = (±L/2, ±L/2). Therefore, the saturated absolute dispersion lies in the range Depending on their initial position, the particles might saturate the domain differently in the spatial directions. Consider, for instance, that the particles are initially released next to the right wall of the large square (L = 1), so (x 0 , y 0 ) ≈ (0.45, 0). Figure 6 presents the particle trajectories at two times before filling the whole domain, in a single numerical run. From panel b, it is noted that the particles overfill the square in the y-direction first, because the horizontal walls are at a shorter distance (L/2) than the opposite vertical wall (located at ∼ L).
To analyse this anisotropy, let us split the saturated absolute dispersion as ⟨a 2 s ⟩ = ⟨a 2 s,x ⟩ + ⟨a 2 s,y ⟩. The saturated components can be calculated as: The saturation value of the total absolute dispersion and its components are tested for the previous example. Figure 7a shows the absolute dispersion curve, which reaches the saturated value ⟨a 2 s ⟩ ≃ 0.3691 (horizontal line) predicted with equation (7). Figure 7b presents the numerically calculated components (9) and (10). They eventually separate from each other, as expected, and reach the saturated values computed with the corresponding formulas.

A. Saturated relative dispersion
The relative dispersion is defined as the mean-squared separation of particle pairs: (11) where N pair is the number of pairs initially separated a distance r 0 and x p,q , y p,q are the positions of particles p and q. 12 Relative dispersion is the second moment of the separation distributions and provides information on the evolution of a cloud of particles regardless of their initial positions. For unbounded, 2D turbulence there are three relative dispersion regimes that depend on the initial pair separation r 0 and the forcing scale r f . These regimes are 6 The first one is the Kraichnan-Lin regime, which occurs for small separations in the range of the direct cascade of enstrophy. 6,31 The spreading of pairs is dominated by turbulent structures with scales larger than the initial separation (nonlocal dispersion), and therefore particles separate exponentially. 32 In the intermediate Richardson-Obukhov regime, the particle separations are in the range of the inverse energy cascade range, so the particles spread under the influence of the energy-containing structures (local dispersion). 6,33 In the third regime, the particles disperse randomly because their separation is larger than the energy-containing eddies. This is the standard diffusion regime where the dispersion is a linear function of time.
Analogously to our previous discussion on absolute dispersion in a confined box, the pair dispersion is saturated as the particles fill the domain. The saturated relative dispersion ⟨r 2 s ⟩ is calculated by averaging the squared separation r 2 = (x p −x q ) 2 +(y p −y q ) 2 between particles p and q uniformly distributed within a square domain: This quantity (1/3 of the area of the square domain) is two times the saturated absolute dispersion (1/6), as expected from the formal definitions of both statistics. 6 Note that, in contrast with absolute dispersion, the saturation value does not depend on the initial location of the particles.
Formula (13) is tested in our numerical simulations. Pairs are formed by using the square array of M × M particles (with M = 32) initially placed at the center of the domain (see Figure 3). We consider three different, dimensionless initial separations q 0 = r 0 /r f = ARTICLE scitation.org/journal/adv i∆/r f , with i = 1, 7, 30 (recall that ∆ = 0.0021 is the minimum separation and r f is the forcing scale). The number of pairs for each q 0 is given by N pair = 2M(M − i), which gives Npair = 1984, 1600 and 128. In all cases there are enough pairs for statistical purposes. Then, the time evolution of the relative dispersion and the kurtosis (the fourth moment of pair distributions) are evaluated. As we did for the absolute dispersion experiments, pair data are obtained from 20 simulations, so the statistics are calculated with 20N pair pairs. All the calculations are repeated for each square domain. Figure 8 shows plots of relative dispersion for sets of pairs with different initial separations q 0 (panels a to c) in the numerical simulations with large, medium and small square domains. Note that the square of the forcing scale (r 2 f ) and the saturation values (L 2 /3) for the three boxes are indicated in all panels. At early times, the three curves of each domain are indistinguishable, which reflects that the dispersion is not affected by the boundaries at this stage. At later times the growth of relative dispersion decreases until reaching the saturation value for each box. The auxiliary curves in the central panel suggest that the exponential and the t 3 regimes take place when the initial separations are smaller than the forcing scale q 0 < 1 (panels a and b). In these two cases the box is saturated during or slightly before the standard dispersion regimes take place. For q 0 > 1 (panel c) the theoretical regimes are not clearly observed, except perhaps the linear standard diffusion regime for the large box. In this case, the relative dispersion is rapidly saturated due to the larger initial separations.

B. Kurtosis of saturation
Another pair statistics is the normalized kurtosis defined by This expression represents the fourth order moment of the pair separation distributions, normalized with the relative dispersion. For the Kraichnan-Lin regime in unbounded 2D turbulence this quantity grows exponentially in time, just as relative dispersion. The kurtosis approaches asymptotic values for the Richardson-Obukhov (K = 5.6) and standard diffusion (K = 2) regimes. 34 For a finite domain, we calculate the limit value of the kurtosis as the particles fill the container. First, we need to evaluate the following four-fold integral: where the integrand represents the fourth power of the separation r between particles p and q. Thus, the normalized kurtosis of saturation is where the saturation value of the relative dispersion (13) was used. This is an exact result that, interestingly, is slightly smaller than the kurtosis of the standard dispersion regime. Figure 9 presents the temporal evolution of the kurtosis calculated in the numerical simulations for different initial separations in each square domain. In all cases, the kurtosis grows at early times (starting at 1), reaches a maximum and then decays at later stages. This behavior has been observed in different studies. 16,33,35,36 The initial growth is more pronounced for q 0 < 1 (panels a and b) because in these cases the exponential Kraichnan-Lin regime takes place (as shown in Figure 8). As the particles separate from each other, the kurtosis is reduced until reaching the saturation value 1.7, predicted in (16). Since the kurtosis is normalized, this result is the same for any square container with arbitrary size. The saturation kurtosis is smaller than the asymptotic kurtosis in the standard dispersion regime for an unbounded domain; therefore, Ks can be regarded as a quantitative measure of the influence of the confinement on the particle statistics. Note that the kurtosis does not reach the asymptotic value of the Richardson-Obukhov regime, expected at intermediate times. This has been a matter of debate in recent dispersion studies. 16,36,37 The present analysis, however, is focused on the modification of the particle statistics and the saturation values due to the flow confinement, rather than the functional form of the time-dependent curves.

IV. SATURATION VALUES IN DOMAINS WITH DIFFERENT GEOMETRIES
In last two sections, saturation values of one and two-particle statistics were analytically obtained by evaluating simple integrals in a square domain. Now we calculate saturation values for different geometries. We consider two sets of shapes: regular (equilateral triangle, square, and circle) and irregular (isosceles triangle, rectangle, and ellipse). Of course, regular domains are special cases of the irregular shapes. The definition of each geometry is shown in the first column of Table II. The geometrical center of each domain is located at the origin of a Cartesian coordinate system.

A. One-particle statistics
The saturated absolute dispersion is calculated first. Consider a set of particles that are homogeneously distributed within a domain with area A. The probability density function (PDF) of particle positions is P(x, y) = 1/A. Thus, the absolute dispersion of saturation for any geometry is defined by where a = (x − x 0 ) 2 + (y − y 0 ) 2 is the distance between the particle position (x, y) and a particle source at (x 0 , y 0 ) (which can be either at the center of the domain or any other arbitrary location). We solve integral (17) for the different geometries. The resulting ⟨a 2 s ⟩ for particles that are initially released at the origin are summarized in the second column of Table II. The saturation values for regular shapes depend on the characteristic length of the domain, e.g. the side length l of the triangle or the radius R of the circle. For irregular domains, the saturation values depend on two characteristic lengths of the geometry, for instance, the semiaxes a, b of an ellipse. Besides, the absolute dispersion of saturation is increased when the original position of the particles is off the origin (third column in Table II), as we saw for the square domain.

B. Two-particle statistics
Now we calculate the two-particle statistics of saturation for the different geometries. Let us consider pairs of particles p and q with uniform distribution over the whole domain, such that the position PDFs are P(x p , y p ) = P(x q , y q ) = 1/A. The joint PDF for the position of a particle pair is The statistical moments of pair separations r = (x p − x q ) 2 + (y p − y q ) 2 are obtained by evaluating the integral: ⟨r n s ⟩ = ∬ A r n Pp,qdA p dA q = 1 A 2 ⨌ A r n dx p dy p dx q dy q .
Using n = 2 the relative dispersion of saturation ⟨r 2 s ⟩ is calculated. With n = 4 we calculate the fourth moment, which is used to obtain the kurtosis of saturation Ks = ⟨r 4 s ⟩ ⟨r 2 s ⟩ 2 .
The values of saturation for both statistics are presented in Table II (fourth and fifth Table II, note a decreasing transition of the kurtosis in regular shape domains: from Ks = 1.8 for the equilateral triangle (a three-sided polygon), to Ks = 1.7 for the square, and to Ks = 1.666. . . for the circle. This suggests that Ks decreases as the number of sides of the polygon is increased, and asymptotically tends to the value obtained for the circle. To test this point, we calculated numerically the kurtosis for a pentagon and a hexagon with random points filling those specific shapes. The results were Ks ≃ 1.678 for a pentagon and Ks ≃ 1.672 for a hexagon, which confirms that the normalized kurtosis for regular shapes lies in the range ARTICLE scitation.org/journal/adv Additionally, we determined the upper and lower limits of Ks for irregular geometries through the following procedure. Consider the kurtosis of saturation for the rectangle with height H and width W (Table II,  (1 + α 2 ) 2 .
By the symmetry of the rectangle, Ks(α) = Ks(α −1 ). Setting the derivative of this expression equal to zero, it is found that the minimum value of the kurtosis is min Ks = 1.7 for α = 1, which corresponds to the square domain. In addition, the kurtosis is maximum if one side differs considerably from the other one, i.e. when H ≫ W (α → 0), or H ≪ W (α → ∞): These limits are verified in Figure 10, which presents the kurtosis in the rectangle as a function of α (solid line). Thus, the normalized fourth moment lies in the range Applying the same procedure to the ellipse with semi-axes a and b, we substitute a = αb into the expression of Ks (Table II,  Ks(α) = 5 6 3 + 2α 2 + 3α 4 and find the extrema of (25) (Figure 10, dashed-dotted line). The results show that the minimum value is min Ks = 5/3 for α = 1, i.e. the kurtosis of a circle. The maximum value (max Ks = 2.5) takes place when the length of a semi-axis is greater than the other, i.e. when a ≫ b (a ≪ b) as α → 0 (α → ∞). Therefore, the kurtosis is bounded in the range 1.666 . . . ≤ Ks ≤ 2.5 for an ellipse.
Finally, the procedure is repeated for the kurtosis of the isosceles triangle of height h and base B (Table II, The minimum value of (27) is min Ks = 1.8, which is the kurtosis for α = √ 3 2, i.e. that of an equilateral triangle with h = ( √ 3 2)B ( Figure 10, dashed line). The maximum kurtosis (max Ks = 2.7) is found when h ≫ B (h ≪ B) as α → 0 (α → ∞). Hence, the normalized kurtosis ranges between 1.8 ≤ Ks ≤ 2.7 for an isosceles triangle.
Let us remark two essential facts from the previous results. First, for irregular domains, the normalized kurtosis approaches a minimum when the two characteristic lengths are such that the geometry approaches a regular shape. Second, the maximum kurtosis is obtained when both characteristic lengths are significantly different, that is, if the size of one length greatly exceeds the other.

V. CONCLUSIONS
The Lagrangian statistics of bounded 2D turbulence were described through numerical simulations. To determine the effects of the flow confinement, we performed several experiments of particle dispersion with square domains of different sizes, in which the characteristics of the dispersive flow were similar. The numerical results showed how the time-dependent curves are modified as the particles fill the domain. At sufficiently long times, the statistical measures reached a constant (saturation) value.
The saturation values for different flow geometries were calculated analytically, and the resulting formulae are summarized in Table II. To the best of our knowledge, these expressions are entirely new. Their relevance stems from the fact that allows predicting the long-term behaviour of dispersion statistics in different closed domains. Also, the present approach provides a theoretical basis for studying more complex, non-symmetrical geometries. The absolute dispersion of saturation depends on the characteristic scales of the domain geometry and on the initial location of the particles, whilst the relative dispersion of saturation only depends on the characteristic scales of the container. The normalized kurtosis of saturation has an exact value for any size of the regular domains. In contrast, the kurtosis depends on the characteristic scales of the irregular containers.
These results can be used to obtain useful information on particle dispersion in a confined space. For instance, we can define a saturation time ts as the time at which the absolute dispersion reaches the saturation value, ⟨a 2 (ts)⟩ = ⟨a 2 s ⟩. The calculation of ts may have practical applications for studying dispersion from an arbitrary site within closed basins (lakes or inland seas). For instance, it provides an estimation of the time at which a pollutant discharged from a point source covers the whole basin. An example based on the present simulations is provided and thoroughly discussed in Flores Ramírez. 38 Another useful result is the modification of the time-dependent statistics as the particles fill the domain and reach a limit value, in comparison with the traditional results for unbounded flows. This is particularly interesting when comparing the saturation value of the kurtosis Ks with the asymptotic value of the random-walk regime K → 2. In our results, Ks < 2 for regular geometries, while Ks > 2 for some elongated shapes of irregular domains (see Figure 10). Previous studies have measured two-particle statistics with surface drifters in the ocean, with the intention to evaluate the presence of the asymptotic limit K = 2. 16,36 Such studies might find important to consider the effects of solid boundaries, as we have shown how they affect the Lagrangian statistics.