On local isotropy and scale dependence of pair dispersion in turbulent canopy flows

Abstract Canopy flows in the atmospheric surface layer play important economic and ecological roles, governing the dispersion of passive scalars in the environment. The interaction of high-velocity fluid and large-scale surface-mounted obstacles in canopy flows produces drag and causes intense, inhomogeneous and anisotropic turbulence. In this work, we focus on the turbulent dispersion of passive scalars by studying the ‘pair dispersion’ – a statistical measure of relative motion between particles. We analyse the results of a three-dimensional particle tracking velocimetry experiment in a wind-tunnel canopy flow, focusing on small scales. We confirm the existence of local isotropy of pair dispersion at scales smaller than a characteristic shear length scale $L_\varGamma =(\epsilon /\varGamma ^3)^{1/2}$, where $\epsilon$ and $\varGamma$ are the mean dissipation rate and shear rate, respectively. Furthermore, we show that pair dispersion in this locally isotropic regime is a scale-dependent super-diffusive process, similar to what occurs in homogeneous isotropic turbulent flows. In addition, we measure the pair relative velocity correlation function, showing that its de-correlation occurs in the locally isotropic regime, and discuss the implications of this observation for modelling pair dispersion. Thus, our study extends the fundamental understanding of turbulent pair dispersion to the anisotropic inhomogeneous turbulent canopy flow, bringing valuable information for modelling scalar dispersion in the atmospheric surface layer.


Introduction
Transport, dispersion, and mixing are essential processes that determine a myriad of physical phenomena with crucial importance to the environment.These processes are driven by the random advection induced by turbulent flows.In particular, numerous such processes occur in the so-called canopy flows, namely in the lower part of the atmospheric surface layer.Examples include the dispersion of viral particles and fungal spores, the ventilation of urban air pollution, and the facilitation of vegetation transpiration by inducing humidity, CO 2 and heat fluxes.Our focus in this study is on pair dispersion in turbulent canopy flows, which involves studying the relative motion between Lagrangian fluid particles.As was shown by Batchelor [4], the statistics of pair dispersion can be used to determine the variance of concentration fluctuations of advected passive scalars.
The study of canopy flows has garnered significant research attention over the past five decades, leading to a plethora of insightful discoveries.This field has witnessed substantial efforts, as evidenced by numerous comprehensive reviews that have effectively summarized the key findings [5,16,17,29,53].
In neutrally stable conditions, canopy flows exhibit a distinctive feature characterized by the direct interaction between fluid and large-scale wall-mounted obstacles within a boundary layer.This interaction generates a drag force that decelerates the flow within the canopy layer (z < H, where H represents the canopy height).At the upper boundary of the canopy (z = H), the drag discontinuity results in pronounced mean shear (∂U/∂z), giving rise to coherent Kelvin-Helmholtz structures akin to a mixinglayer analogy [52].Turbulence statistics within the canopy layer are consequently both inhomogeneous and anisotropic.Eulerian single-point velocity probability density functions (PDFs) typically exhibit skewness, indicating intermittency and relatively high kurtosis values (Brunet, 1994).Interestingly, Lagrangian velocity increments have been observed to follow a Gaussian distribution [68].Additionally, turbulence in canopy flows is typically generated through two interdependent mechanisms: (I) shear production near z ≈ H, and (II) production of smaller-scale turbulence within the wakes of the obstacles within the canopy (z < H) [29].
A common framework for modeling dispersion in canopy flows is through Lagrangian stochastic models [?, 1-3, 21, 24, 27, 28, 31, 45, 51, 54].In such models, Lagrangian fluid particles are advanced through the flow field using random walks that simulate the turbulent flow; this allows estimating statistics of the concentration field as a result of a certain distribution of scalar sources.Although these types of models allow us to bypass such difficulties as the hindering of parameterizations due to the existence of multiple scales or the inaccuracy of Taylor's frozen-turbulence hypothesis [51,53], there are still numerous open issues in their development and formulation.For the case of single particle motion, these issues include the non-uniqueness of first-order Markov random walk models in three dimensions [75], the effects of coherent structures [30,52], non-Gaussian velocity PDFs [48], the parallel contributions of wake and shear production [68], or the mechanical diffusion [38].To the best of our knowledge, previous studies have only considered the motion of single particles, and there are no previous studies that deal with the relative motion between particles in canopy flows.Notably, a full description of the dispersion of a group of particles demands that the relative motion between any of their possible combinations be resolved [4].
Pair dispersion deals with the relative motion of pairs of particles, thus making a step forward in this respect.
In this work, we present an analysis of pair dispersion in a canopy flow using results from a wind tunnel experiment.This flow mimics the neutrally buoyant atmospheric surface layer in an environmental windtunnel setup, as we reported in Ref. [66].We focused in the past on single particle statistics of this canopy flow, revealed short decorrelation Lagrangian timescales [68], and characterized the unique features of Lagrangian intermittency in this flow [67].Our results deal with relatively small scales and are relevant for describing dispersion in between roughness elements inside the canopy.These findings are useful for instance for modeling dispersion between buildings in urban areas.
Background on pair dispersion and details on the experiment are given next.Following that, our analysis is divided into two main parts -first we quantify the effect of the mean shear and anisotropy in our canopy model through the pair dispersion tensor.Following that, we study the different regimes of pair dispersion that are observed in our measurement, including the ballistic and the inertial regime with a scale-dependent diffusivity; we suggest how these results could be useful for estimating the concentration variance of released passive scalars in the flow.Thus, our results could be useful for constructing and validating turbulent dispersion models in the atmospheric surface layer.

Background on pair dispersion in isotropic and anisotropic turbulence
Turbulent pair dispersion was introduced by Richardson [55] and since then it has been studied extensively, e.g. as reviewed in Refs.[26,56,57].Pair dispersion describes the statistics of vector r connecting two, initially close, Lagrangian fluid particles: where x (i) is the 3D coordinate of a particle i (bold symbols denote vectors).Leaning on principles of homogeneous isotropic turbulence (HIT), the dynamics of r can be divided into three distinct regimes.(I) At short times, there is a ballistic regime [4], during which particles' relative motion is co-linear.(II) At larger time scales, Richardson [55] suggested that statistics of r can be modeled as a diffusive process with a scale-dependent diffusivity parameter, K: where ⟨•⟩ is an ensemble average of many particle pairs.The scaling of K, commonly termed "the fourthirds law", leads to a super-diffusive growth of the separation distance with r 2 ∼ τ 3 [37].(III) at later times when the two particles have separated farther away than the integral turbulence scale, i.e. for r ≫ L, the distance r is expected to be diffusive with a constant, or Taylor, diffusivity [73]: where T L is the Lagrangian timescale and ũ is the root-mean-square (RMS) of velocity fluctuations.
Therefore, considering an ensemble of pairs of particles, initially separated by r 0 and moving in a HIT flow, the second order moment of the distance, the variance of the change in separation distance is summarized as follows [37,56]: where C 2 ≈ 2.1 [70] is the so-called Kolmogorov constant, ϵ is the mean rate of the turbulent kinetic energy dissipation, g ≈ 0.5 [41] is a so-called universal Richardson constant, and the Batchelor timescale is defined as: The theory of pair dispersion has been critically revised and extended in recent years.In particular, the three regimes have been tested in several experiments and numerical simulations for various initial distances, r 0 /η (where η = (ν 3 /ϵ) 1/4 is the Kolmogorov dissipation length scale).Verification tests of Eq. ( 4) in Refs.[6,8,11,15,41,43,77] have shown that it is difficult to observe the transition from ballistic to the super-diffusive regime, partly because it requires long-duration tracking, i.e. for τ ≫ τ b , and high Reynolds number flows [6,15,42].Specifically, the super-diffusive regime is asymptotic, and so, finite Reynolds number effects cause the evolution of statistics of r to depend strongly on the initial separation, r 0 /η [8].Another issue is the strong intermittency of pair dispersion [9,10,58,65,74], presumably due to Lagrangian statistics associated with long-time correlations [12].In particular, recent observations have shown that the scaling of r 2 depends on the initial conditions, namely on the initial separation [25,72] and the initial separation velocity [65], while it is also known that particles may retain their separation distance for very long times [58].Thus, it was suggested that Richardson's diffusive approach relies on an assumption of short Lagrangian correlation times [10], which could explain the inconsistency with the recent observations.Several models using different forms of the so-called persistent ballistic random walks have been proposed as alternatives [14,50,69,74].Bitane et al. [10] proposed an alternative to τ b which renders the process self-similar with respect to r 0 ; this framework is discussed in more detail below.
While Eq. ( 4) and the studies mentioned above pertain to HIT flows, pair dispersion in inhomogeneous or anisotropic flows has received much less attention.In inhomogeneous and anisotropic flows, statistics vary for each component of r, and their evolution may be co-dependent on each other.In particular, in anisotropic flows it is necessary to introduce a pair-dispersion tensor with components [4]: notably, tr(∆ ij ) = (r − r 0 ) 2 , and in HIT ∆ ij is diagonal.Furthermore, pair dispersion in anisotropic flows can depend on the initial orientation of r 0 with respect to the production mechanisms of turbulent kinetic energy; for shear flows, such as the canopy flow considered here, this would be the mean shear direction.
Shear effects on pair dispersion were studied previously in homogeneous shear flows through computational simulations.Shen and Yeung [61] established that the presence of a mean shear has a strong effect on pair dispersion: it makes the evolution of ∆ ij anisotropic and with non-zero non-diagonal components, and it can render apparent super-diffusivity at long times (∆ xx ∼ τ 3 ) irrespective of the argument of scale-dependence.More recently, Pitton et al. [44] and Polanco et al. [47] studied pair dispersion in turbulent channel flows and highlighted the combined effects of both anisotropy and inhomogeneity with respect to the distance from the wall.Their analysis clearly showed that where the mean shear is strongest the anisotropy of pair dispersion is most dominant.In a DNS of a square duct flow, Sharma and Phares [60] reported that particles remained trapped in secondary flow regions which affected the power law behavior of r(τ ).Celani et al. [19] have proposed that mean shear and turbulent fluctuations act on pair dispersion on different time scales, which leads to two distinct regimes of either turbulent fluctuations dominance or mean shear dominance.A crossover between these two regimes was predicted to occur at a critical timescale that is proportional to the mean shear τ c ∼ dU dz −1 .Recent studies further examined buoyancy-driven flows [34,39,59] and pair dispersion of inertial particles [44] or in other complex flows, such as in the stroke of a swimming jellyfish [32].
Despite the extensive body of work, the above discussion establishes the existence of a wide gap in the understanding of turbulent pair dispersion.Indeed, even for ideal cases, such as HIT or homogeneously sheared turbulence, our understanding is lacking and the introduction of flow inhomogeneity brings on another dimension of complexity.In this study, we bring the first empirical results on pair dispersion from an inhomogeneous and anisotropic canopy flow in order to shed some light on the prevailing issues.

A 3D particle-tracking wind-tunnel experiment
We study pair dispersion in a canopy flow using the results of a 3D particle tracking velocimetry experiment (3D-PTV [23]).In the experiment, we obtained flow tracer trajectories in a wind-tunnel canopy flow model using an extended, real-time image processing system [66].The wind tunnel laboratory, situated at the Israel Institute for Biological Research (Ness Ziona, Israel), features a 14 meters long open circuit suction wind tunnel with a 2×2 squared meters cross-sectional area, that is fit for conducting experiments that mimic turbulent flows in the atmospheric surface layer [13].The canopy flow was modeled by placing flat rectangular plates along the entire bottom wall of the test section.Our study used an inhomogeneous canopy layer, constructed with two types of flat plates with heights either H or 1 2 H, and width 1 2 H, where H = 100mm; the two types of plates were positioned in consecutive rows and at a staggered construction, see Fig 1 .The canopy frontal area density, defined as Λ f = A f /A T , (A f being the element frontal area, and A T the lot area of the canopy) is 1 2 , which positions our canopy right between the "dense" and "sparse" categories [17].Data was gathered at two levels of the free stream velocity, U ∞ = 2.5 and 4m s −1 being the free stream mean velocity measured with a sonic anemometer at the center of the wind tunnels' cross-section.These values correspond to Reynolds numbers Re ∞ ≡ U∞H ν ≈ 16 × 10 3 and 26 × 10 3 , where ν is the kinematic viscosity of the air.
Our previous work with this canopy construction, based on particle image velocimetry (PIV) measurements, showed that the double-averaged velocity profile in our canopy layer had a smaller shear length at the top of the canopy as compared to homogeneous canopy layers at similar densities and that, unlike homogeneous canopies, it had a second inflection point right above the lower roughness elements [63].
Furthermore, as commonly observed in canopy flows, sweep events were observed inside the canopy layer (z < 1.5H), albeit with a second peak of increased sweep contribution above the lower set of elements [62].
In addition, our previous analysis of 3D-PTV measurements showed that the canopy is characterized by an unexpectedly short Lagrangian de-correlation timescale, leading to probability density functions of velocity increments that are of a Gaussian shape [68] and complemented by intermittency, as observed in single particle Lagrangian statistics [67].
In the experiment, we tracked the positions of flow tracers using the 3D-PTV method through an extended 3D-PTV application.The details of the methods are given in Ref. [66], so only brief information shall be repeated here for completeness.The flow was seeded with hollow glass spheres with an estimated Stokes number of St ≈ 0.05, indicating good tracking properties with only a minor filtration of high-frequency contents of the turbulent dynamics that cannot be ruled out.The tracer particles were illuminated with a 10W continuous wave laser beam as shown in Fig 1a .Using a set of four 4-megapixel cameras, we obtained 3D particle trajectories at a resolution of 50µm per pixel, inside volumes of measurement were roughly 8 × 5 × 4mm 3 in size.We implemented camera calibration, stereo matching, and tracking [23], to reconstruct the tracer particle's trajectories by using the openPTV [40] open-source software, integrated to operate with the real-time image analysis extension.The data analysis was performed using the Flowtracks package [36].In this work, we use the frame of reference that is commonly used in the canopy flow literature: x is aligned with the streamwise direction longitudinally within the wind tunnel, y is the cross-stream horizontal direction, and z points vertically up against gravity and away from the bottom wall.

Data analysis: a quasi-homogeneous approach
We recorded trajectories in the height range 0.5H < z < 1.5H and across a single representative canopy unit cell, by scanning the full volume through 20 sub-volumes, as depicted by the green shaded region in Fig. 1a.The sub-volumes were defined at 4 horizontal locations, and at 5 different heights above the wind-tunnel bottom wall.Sub-volume horizontal positions are labeled alphabetically a, b, c and d as shown in Fig. 1(b).At each position, a − d, we used 5 vertical slabs of thickness 0.2H to define the subvolumes; the vertical slab position is labeled numerically 1-5 and corresponds to heights z/H ∈ {0.5-0.7,0.7-0.9,0.9-1.1,1.1-1.3,1.3-1.5},respectively.The volume scanning approach allowed us to measure a full canopy unit cell, however, it limited the particle trajectory lengths to be within each of the individual sub-volumes.
In the frame of our analysis, we present statistical properties calculated over ensembles of trajectories.Thus, formally, the mean value of any Lagrangian quantity A is estimated using the average - where V is a tag pertaining to different particle ensembles.In our recent work [68], we explored singleparticle statistics using the same dataset.There, the sub-volume averaged quantities approach was proposed, in which Lagrangian statistics were presented for ensembles of particles found at each sub-volume, essentially resulting in spatially averaged statistics across each small sub-volume.It was furthermore shown that this quasi-homogeneous approach is justified in our flow because random turbulent forcing out-weighted combined effects of flow inhomogeneity terms in the stochastic particle equation of motion.
Specifically, this was achieved through the condition: 1 2 where C 0 is the so-called Lagrangian structure function coefficient, R ij = u ′ i u ′ j is the Reynolds stress tensor, and ϕ i /g is a vector representing a sum of drift terms of the Lagrangian velocity PDF that allows accounting for flow inhomogeneity [75].In the present work, we capitalize on the condition (8), and present quasi-homogeneous statistics for pair dispersion as well.This point is important since it validates the effective scaling laws we obtain below.
In addition to that, in Ref. [68] turbulent parameters of our flow were calculated where they were discussed and analyzed in depth.In particular, we calculated ϵ, the turbulence dissipation length and time scales η, τ η , the Lagrangian decorrelation timescale T i and the Taylor Reynolds number Re λ .The values of these parameters are tabulated for each sub-volume in Appendix A, and they shall be used in the analysis below.

Local isotropy of pair dispersion 3.1 Scaling analysis
As discussed in Sec. 1, numerical simulations [44,46,47,60,61] confirmed theoretical predictions [37] that the mean shear, inhomogeneity, and/or anisotropy in the flow cause anisotropic pair dispersion.It is based on the fact that the fundamental mechanism driving pair separation is the relative velocity between Lagrangian particles.Because turbulent flow statistics are scale-dependent, it is expected that also the "degree" of anisotropy in pair dispersion will depend on the separation distance r.Turbulent velocity fluctuations are characterized by a tendency to recover isotropy at small scales due to the isotropy of the equation of motion.The return to isotropy at smaller scales was also observed in canopy flows [71].This suggests that at sufficiently small values of r, and as long as the Reynolds number is sufficiently high, pair dispersion will be, at least approximately, isotropic, as predicted theoretically in Ref. [19].Thus, to determine whether local isotropy or anisotropy is expected to occur in our pair dispersion measurement, we first conduct an analysis of the relevant scales of the flow.
Anisotropic turbulent fluctuations are introduced to the flow by the boundary conditions.In the canopy flow case, the flow impinges on the roughness obstacles that exert drag and retard the flow.Thus, the production of turbulence in canopies is commonly decomposed into two main contributions [29]: I) shear production at the top of the canopy layer, and II) production at smaller scales in the wakes of the obstacles.The third possible component, the production by wall friction is usually negligible compared to form drag due to the low velocity near the ground for sufficiently dense canopies.Because shear production relates to the horizontally averaged mean shear, it induces an intrinsic anisotropy in the structure of turbulence in canopy flow.On the other hand, wake production occurs due to local variations in the rates of strain, associated with the wakes and boundary layers of individual roughness element [29].
Since the local orientation of most straining directions changes in space, we expect that the anisotropy in pair dispersion statistics will mostly reflect shear production, at least to a first approximation.This seems to be particularly true for Lagrangian statistics such as in pair dispersion since the particles essentially sample different flow regions.
The appropriate length scale that characterizes the production of turbulence by shear is where Γ = dU dz is the mean velocity gradient.Indeed, analyzing the scale-by-scale turbulent kinetic energy budget in turbulent shear flows, Casciola et al. [18] showed that the production by shear continues down to L Γ , whereas for scales r < L Γ the TKE transport becomes more dominant.For our analysis, we parameterize the shear length using a global, canopy-wide measure, of the mean velocity gradient - 2(a) shows that L Γ /η increased from roughly 100 inside the canopy and up to roughly 150 above it.
Another critical dimension of pair dispersion is time.Celani et al. [19] predicted that for particles with sufficiently small initial separations, pair dispersion is initially isotropic while it becomes anisotropic once the particles have grown sufficiently far apart for the anisotropic turbulent scales to be prominent.
They further noted that the critical timescale for this transition is proportional to the inverse mean shear rate, τ c ∼ Γ −1 .Figure 2(b) shows the Lagrangian velocity decorrelation timescale for the streamwise velocity component normalized by Γ.Throughout our measurement region, we observe that T x • Γ < 1.
Due to the finite measurement region, our work focuses on timescales up to τ ≈ T x .In addition, the typical separations between particles that we consider in our work reach up to roughly r = 100η, namely mostly values with r < L Γ .Thus, the dimensions of our measurement region and the timescales on which we focus are smaller than the scales imposed by the mean shear.For this reason, we expect to see only

Observation of weak anisotropy
The various components of the pair dispersion tensor, ∆ ij , across the whole region of measurements are shown in Fig. 3 for two representative initial separation values inside the inertial range.In all cases, the diagonal components steadily increase from zero as time increases, which indicates increasing separations in all orthogonal directions as expected.The off-diagonal components, on the other hand, generally do not increase significantly in the available time frame of the measurements, with a trend for larger scatter for the larger initial separation cases and at longer times.This occurs consistently across the entire measurement sub-volumes.
Since the non-diagonal components represent mixed averages of the separation vector's components, their relatively low values indicate a weak correlation among the separation in different directions.This behavior is unlike what is expected in a mean shear-driven separation scenario as, for example, in a homogeneously sheared turbulence the separation velocity in the streamwise direction grows as the separation along the mean velocity gradient direction increases.Thus, the behavior observed for ∆ ij is consistent with the notion presented in Sec.3.1 that the mean shear is not expected to significantly affect pair dispersion in our measurements.
The three diagonal components of the pair dispersion tensor grow with time, and slightly different values are seen for the various components.To highlight these differences we introduce the following tensor where δ ij is the identity matrix.The tensor I ij shows statistics of the components of ∆r with respect to its norm, essentially highlighting anisotropy.Let us note that tr(I ij ) = i I ii = 0, and that dispersion in the direction of components with I ii > 0 is faster than in the isotropic case (i.e., than 1 3 (r − r 0 ) 2 ), while components with I ii < 0 the dispersion is slower.In a limiting case in which pairs separate only in one direction, say x, then I xx = 2 3 while I yy = I zz = − 1 3 .Histograms of the diagonal components, I ii (τ ) are shown in Fig. 4, focusing on two representative initial separations and three representative heights: inside the canopy 0.5 < z/H < 0.7, at its top 0.9 < z/H < 1.1, and above it 1.3 < z/H < 1.5.The histograms count data points from all horizontal The disparity between the separation components decreases with height, being the largest inside the canopy and weakest above it.Furthermore, the separation was the slowest in the vertical direction and fastest in the spanwise direction, whereas the streamwise separation is approximately at isotropic values (I xx ≈ 0).The fact that the separation is fastest in the spanwise direction and inside the canopy, suggests that the observed anisotropy could be attributed to the channeling of the flow in between the roughness obstacles, which were in a staggered configuration, and not necessarily due to the effects of the Figure 4: Diagonal terms of the pair dispersion tensor normalized by its trace minus one-third.Data is shown for sub-volumes at three heights and for two levels of r 0 .mean vertical shear.This would highlight the effects of dispersive fluxes between the canopy obstacles in driving horizontal dispersion.And still, we note that the magnitudes of I ii values reach only up to approximately 0.1; this corresponds to a weak anisotropy with 20% of disparity between the separation in the vertical and spanwise components.In particular, this degree of anisotropy is much smaller than what was observed in previous works concerning flows with mean shear [44,47,61].
The disparity between the components is also seen to increase with the initial separation, r 0 .This is consistent with the picture of return to isotropy, as anisotropy becomes more prominent at larger scales [35].
The increasing of anisotropy with scale brings up questions regarding the development of anisotropy in pair dispersion with time.As particles separate, the growth of typical r values suggests that pairs are exposed to the more anisotropic turbulent scales as time increases.The development of anisotropy could be examined by adopting a framework analogous to the one introduced by Stiperski et al. [71], which used a projection of the invariants of the normalized Reynolds stress tensor on a two-dimensional plane, allowing them to investigate trajectories of the return to isotropy in various canopy flows.In an analogy to that, we examine here the projection of the eigenvalues of I ij using the same projection: where λ 1 , λ 2 , and λ 3 , are the smallest intermediate and largest eigenvalues of I ij respectively.Equation 11maps the eigenvalues to a planar equilateral triangle.As seen in Fig. 5, the triangle's nodes correspond to cases of fully isotropic, 2-component axisymmetric, or one-component dispersion.Plotting the trajectories that correspond to the measured I ij tensor allows to probe the topology of pair dispersion as it varies in time.For almost all of the data shown, the weak anisotropy of pair dispersion is seen by the fact that the trajectories are almost completely confined to the isotropic part of the map.Nevertheless, the weak anisotropy that does exist is evident in the general trend of trajectories to progress along the left flank of the triangular maps, in the direction corresponding to the oblate topology.The trajectories also demonstrate that most of the anisotropy in our pair dispersion measurements is explained by the larger initial separation values and not due to the increase in separation with time.This is in line with the fact that our measurements are confined to time scales smaller than Γ −1 .

Bias due to initial orientation
As shown in Refs.[44,47,61], in flows with mean shear, pair dispersion is affected by the initial orientation of the pairs for short times.Indeed, taking as an example the case of a shear flow where the velocity is u = (Γ z, 0, 0) and considering particles with an initial separation r 0 = (r 0,x , r 0,y , r 0,z ), then the separation vector is r(τ ) = (r 0,x + τ Γ r 0,z , r 0,y , r 0,z ).Therefore, even in the simplest scenario, the presence of mean shear biases the streamwise separation component.The bias depends on the projection of the initial separation vector on the mean shear direction, so for pairs whose r o is aligned with the mean shear direction, the bias is by a factor of Γr 0 .
The relative importance of this bias in a turbulent flow can be estimated by comparing Γr 0 with the spatial fluctuating velocity increments δ r0 u.If we apply the Kolmogorov dimensional scaling [33] to parameterize δ r0 u ∼ (ϵr 0 ) 1/3 , we can construct a dimensionless group S = Γ r0 δr 0 u = Γr 2/3 0 that quantifies the importance of the shear bias for particles with r in the inertial range.When S ≫ 1 we expect that shear will bias pairs separation based on their initial orientation whereas when S ≪ 1 it will not.Also, if an inertial range exists, this bias grows stronger as r 2/3 0 for smaller initial separations, which is in qualitative agreement with the observations of Refs [47,61].Also, in the dissipation range, the velocity field is presumably smooth and the typical size of velocity gradients is 1 τη = ( ϵ ν ) 1/2 (where τ η is the dissipation time scale).Therefore, the dissipation scaling for velocity increments in the dissipation range is δ r0 u ∼ r 0 ϵ ν 1/2 so that the dimensionless parameter becomes S = Γ ϵ ν −1/2 .To summarize, we define and when S ≫ 1 the presence of mean shear is expected to bias pair dispersion, making it faster when the initial orientation of particles is aligned with the shear direction.
To examine the effect of initial orientation in the canopy flow, we estimated the variance δr 2 conditioned on the initial orientation of r 0 .We divided the pairs into sub-samples based on the condition that the angle between r 0 and each of the coordinate axes was lower than 25°.In Fig. 6, probability distributions are shown for the conditional variance normalized by the variance of all the pairs.In this case, data was taken at all available times, all sub-volumes, and with 30η < r 0 < 50η.The distributions show that pairs with initial separation aligned vertically (i.e. with ẑ) typically separated faster than the average, whereas pairs with initial separation aligned with the spanwise (i.e.ŷ) typically separated slower than the average.Considering the previous observation that separations are fastest in the spanwise direction (Fig. 4), the bias observed here for vertically oriented pairs may suggest that the channeling effect varies with height.This is also consistent with the observation of the strongest anisotropy inside the canopy layer.
From Fig. 2 it is seen that for the present case, we have a maximum value of S ≈ 0.38 which suggests that the bias should be quite weak.This is validated in Fig. 6 since it shows that the magnitudes of the pair dispersion bias encountered in our measurements reached up to 30% and was typically around 15%.These values are indeed rather low considering that previous studies [44,47] observed orders of magnitudes of difference between different groups of conditioned particles in regions of the flow with strong shear.Overall, our empirical results agree well with the scaling argument presented above and are consistent with the theory of Celani et.al [19].These results reveal the important role of L Γ and Γ in determining the scales relevant for the development of pair dispersion anisotropy in canopy flows.Indeed, anisotropy was weak at the scales relevant to our work, in addition to a weak bias due to the initial orientation of pairs.These key results highlight local isotropic turbulent fluctuations as the main driver of particle separations in the canopy flow at small scales.

Regimes of pair dispersion and a scale dependent diffusivity
Following the observation in Sec. 3 that pair dispersion is only weakly anisotropic in our measurements we focus next on statistics of the change of the separation distance, r = |r|, with time.In particular, we analyze the scaling of the change in separation distance (r − r 0 ) 2 ∝ τ β [6,8,43,56,65].It is also noted that, formally, due to the inhomogeneity of our canopy flow, β is essentially an "effective scaling".
Therefore, the goal of this section is to compare the empirical data with the different regimes of Eq. ( 4), in order to understand the phenomenology underlying pair dispersion in our canopy flow.
In order to investigate time scaling, it is essential to consider a possible bias due to the finite volume of measurement.In particle tracking experiments there is a concern that finite volumes could lead to a bias in the estimate of Lagrangian scaling exponents since particles leave the volume at different times.
A similar bias when estimating Lagrangian correlation functions was investigated in Ref. [7].Specifically, the extreme cases in which particles separate very fast and leave the volume of measurement at short times could influence the slopes of the curves leading to errors in estimating the scaling exponent β.
To avoid such a bias, the data were sub-sampled for pairs tracked for a fixed amount of time.For our analysis, we consider only pairs that were successfully tracked for durations t tr > 5τ η .Statistics were calculated only for τ ≤ 5τ η (where t tr is the tracking time, the duration of time for a particle to remain in the sub volume).This sub-sampling may cause an underestimation of the rate of separation.However, it is necessary for making reliable estimations of the scaling exponent.In addition to that, to maintain a concise description we present data for only one sub-volume -b3, albeit the results were observed to be robust for particles from the other sub-volumes as well.

The ballistic regime
We first verify that pair dispersion at short times follows a ballistic propagation, as was suggested by Batchelor [4].This regime is characterized by a quadratic growth of ∆r 2 ∝ τ 2 (τ ≪ τ b ), corresponding to the first case of Eq. ( 4).Fig. 7 (a) presents pair dispersion from the canopy flow in normalized form according to Eq. ( 4), where S LL (r 0 ) = (δ r0 u • r0 r0 ) 2 is the directly measured Eulerian structure-function, and the Batchelor timescale τ b .Initially, for τ < 0.1τ b , all the curves with different initial separations collapse on a quadratic dashed line, as Eq. ( 4) suggests.Similar results were previously obtained in homogeneous flows by Refs.[10,43,65], and inhomogeneous flows by Refs.[34,39,47].The collapse of the data at short times using the longitudinal structure function S LL (r 0 ) is important since it confirms that [S LL (r)] 1/2 is the appropriate velocity scale for the separation of particles at a distance r.

Transition from a ballistic to an inertial regime
Next, we examine the termination of the ballistic regime, which is done following the analysis performed recently by Bitane et al. [10] for the case of HIT flow.We denote the separation velocity v || = ∂r ∂t and acceleration a || = ∂v || ∂t .Then, the separation distance is Taylor expanded to the third order, squared, and ensemble-averaged, which gives where subscript 0 denotes initial values at τ = 0.This expression, which extends the Ballistic regime, is expected to hold for short times, and, in particular, according to Bitane et al. [10], it gives the appropriate timescale for the end of the ballistic regime at . Note that the minus sign is used since in 3D turbulence a ||,0 v ||,0 < 0 [41], as was confirmed for our data (not shown for brevity).Therefore, the above Taylor expansion can be rewritten in dimensionless form as In Fig. 7(b) we show the empirical data of ∆r 2 for the different r 0 cases, normalized according to the left-hand side of Eq. ( 15), and compare with the right-hand side shown as a dashed black line.For short times, τ ≪ τ 0 , the figure shows rather good agreement between the empirical data and Eq. ( 15), especially in the downward trend of the lines, which means that the above arguments of Ref. [10] agree with our data.Furthermore, the negative term in Eq. ( 15) is responsible for the deviations from a purely  13).(b) Same data but normalized using the time lag, and the dashed line corresponds to Eq. ( 15).quadratic growth as time progresses.Therefore, defining the end of the ballistic regime as a deviation of ∆r 2 by, say, 10% from v ||,0 τ 2 , we obtain that the ballistic regime ends at τ = 0.2 τ 0 .Indeed, our data from different r 0 agrees with Eq. (15) (dashed line) roughly up to τ ≈ 0.2τ 0 .
Once the ballistic regime is terminated, the separation velocities of particles have deviated significantly from their initial values.This can be shown by examining the autocorrelation of the separation velocity, defined here as which is shown in Fig. 8.The autocorrelation ρ v || is plotted against τ τ0 for the 5 cases of pairs with different r 0 .The graph shows that the correlation of the separation velocity drops to roughly zero once the time lag is τ ∼ τ 0 (with some weak dependence on r 0 ).This confirms that at τ ≳ τ 0 the particles' separation was with velocities that are very different from their initial value.We shall call this regime of the pair dispersion inertial.

Scale-dependence growth rates and approach to Richardson's 4/3 law
Richardson's classical theory [55] describes pair dispersion as an accelerating process in which the rate of separation increases with r.Accordingly, as r increases with time, we should observe that the separation rate grows as well.Nevertheless, as discussed in the introduction, this is often challenging to observe in experiments and simulations due to the need for sufficiently long observation times such that r increases significantly.To confirm whether such conditions occurred in our experiment we show in Fig. 8(b) the normalized separation, r/r 0 , for five groups of pairs with different initial separations.For the group with the smallest initial separation available, r 0 = 10η, the typical separation increased by a factor of two by the end of our measurement range; on the other hand, the increase in separation was much lower for the other groups with larger r 0 .This issue, however, is typical in experiments and simulations [8,25,34,41,72], because the separation timescale, τ 0 ∼ r 2/3 0 , is smaller for smaller r 0 .For this reason, if an explosive pair dispersion regime occurred in our experiment, it should be most evident in statistics of particles with small initial separations.Indeed, in Fig. 7 (a), an increase of the local slope for the r 0 = 10η group can be seen at longer times (τ ≥ τ 0 ), indicating an increase in local scaling exponent.In these cases, the pair dispersion is super-diffusive, i.e. ∆r 2 ∝ τ β with β > 2.
With the observations obtained thus far, let us lay out the reasoning to explain the behavior of pair dispersion in Fig. 7 at longer times, namely after the Ballistic regime is over.First, for all cases of r 0 the separation velocity became decorrelated in our measurement (i.e.Fig. 8(a)), so that a so-called inertial regime was reached.Second, the increase of r relative to its initial value depended on r 0 , as shown in Fig 8(b).Third, in Richardson's theory pair dispersion is described as a diffusive process where the diffusivity depends on r (i.e.K r (r)).In the present case, pairs with small r 0 had multiplied their separation distance within the range of our measurement and therefore had increased their rate of separation, leading to the super-diffusive scaling β > 2. In contrast, for the largest r 0 case there was ⟨r 2 ⟩ 1/2 r0 ≈ 1 for the entire range of our measurements.Therefore, in Richardson's diffusive framework, these pairs had an almost constant diffusivity, so the scaling exponent was β ≈ 1, as expected in a diffusive process with constant diffusivity.We support the above picture by estimating diffusivity directly.Following Batchelor [4], we define In Fig. 9(a) the diffusivity is plotted against the RMS r ≡ r 2 1/2 , where the axes are normalized using the dissipation scales, η and τ η .An arrow on the figure indicates the direction in which the time grows.
For all cases shown, K r is initially negative, which, as was discussed by [49], is a consequence of the negative skewness of the Eulerian velocity differences, aka the 4/5 law [37].With increasing time K r grows.At small times K r grows in line with the Ballistic regime, the end of which, at τ = 0.2τ 0 , is indicated as circles on the figure.Later on, at τ ≈ τ 0 , the separation velocity became decorrelated which is followed by the inertial regime.The time τ = τ 0 is indicated by stars on the figure.The range of our measurement extended to this regime (τ ≳ τ 0 ) for the cases of small r 0 .
In the inertial regime, most pronounced for the small r 0 cases, the curves tend to the right side, meaning that as r increased, the diffusivity increased as well.This is the main mechanism leading to the accelerating pair dispersion picture.Furthermore, the Richardson 4/3 law (as interpreted for the averaged separation by Batchelor [4]), according to eq. ( 2) predicts that K r ∝ r 4/3 .The prediction eq. ( 2) is plotted in the figure as a dashed black line.The data for the two smallest r 0 cases leans towards the dashed line for τ > τ 0 .This observation agrees with previous results that showed a super-diffusive, t 3 scaling range for small initial separations in homogeneous isotropic turbulence [25,72].
The Kolmogorov scaling for the diffusivity in the inertial range is K r ∼ (r 4 ϵ) 1/3 .Therefore, we annotate K r0 ≡ (r 4 0 ϵ) 1/3 .Accordingly, in Fig. 9(b) the same data is plotted where the axes are normalized using r 0 and K r0 .The data from all r 0 cases is seen to collapse on a single curve under this normalization.
The collapse of the data from all r 0 in Fig. 9(b) confirms the existence of a scale-dependent regime in our pair dispersion measurement.This is a key result of this work.Also, as before the data for the smallest r 0 case joins the K r = K r0 (r 4 ϵ) 1/3 line which corresponds to Richardson's 4/3 law, Eq. ( 2).
These observations extend the fundamental phenomenology of turbulent pair dispersion to small-scale separation in highly turbulent anisotropic flows.

Discussion and Conclusions
Our work presents the first empirical investigation of turbulent pair dispersion in a canopy flow.For the analysis, we have used an extensive dataset of Lagrangian trajectories measured in a wind-tunnel experiment using an extended 3D-PTV method [66].
The first part of our analysis focused on the anisotropy of pair dispersion.Our investigation is focused on small scales, being of the order of separation between canopy obstacles.As opposed to the anisotropy of canopy flow in larger scales, our measurements show that in the canopy, at the scales of separation between canopy obstacles, pair dispersion anisotropy is weak.The disparity among the separation components reached only up to approximately 20% of the total separation.This was the result of our measurements being focused on scales smaller than L Γ -the scale that determines the crossover between small-scale local isotropy and large-scale anisotropy [19].In addition, the time range of our measurement was smaller than Γ −1 , so even though particles separate with time and the scales to which they are exposed grow, their separation was still insufficiently large for anisotropy to build up and become prominent.
The weak anisotropy that does exist was observed to grow with time and initial separation.Furthermore, the spanwise dispersion inside the canopy was observed to be the strongest contributor to anisotropy in our measurements, presumably due to flow channeling in between the roughness obstacles.Investigating the pair dispersion tensor in a framework analogous to the return to isotropy [35,71] (Fig. 5), showed that the anisotropy trajectories of pair dispersion progress towards anisotropy on a path that approximates oblate spheroids.Interestingly, in a recent meta-analysis of the Reynolds stress tensors anisotropy trajectories across numerous surface layer flows, a similar trajectory along the oblate topology was observed [71].This similarity might hint that theory developed in the framework of "return to isotropy" could provide helpful information on pair dispersion modeling in anisotropic flows; however, this is left for future work.
An important feature observed in our experiments was the smallness of turbulence scales compared to those of mean shear.Due to this feature, the locally isotropic turbulent fluctuations were sufficiently strong to drive the nearly isotropic pair dispersion.This highlights the importance of L Γ /η and T L Γ to understanding pair dispersion in anisotropic turbulent flows, as having high and low values respectively implies that the isotropic turbulence phenomenology is the appropriate tool for analysis.
The second part of our analysis focused on the next two regimes of pair dispersion using the conceptual framework of local isotropy.At short times, the particles follow ballistic trajectories with separation velocities dictated by the Eulerian structure-function, as observed in homogeneous turbulence [15].At later times, the separation velocity undergoes decorrelation and a so-called inertial regime ensues.The rate of separation was shown to depend on the scale, r.For pairs with small initial separations, the separation rate was seen to scale in accordance to Richardson's 4/3 law (r 0 = 10η), in agreement with previous observations in isotropic turbulence [25,72].These observations confirm that the leading mechanisms that govern pair dispersion at small scales in canopy flows are inherently tied to the internal turbulence regulation where the kinetic energy grows with the scale.
Understanding pair dispersion is important for modeling dispersion in the environment, as it can be used to estimate the variance of passive scalar concentration fields from known sources.In particular, Cohen and Reynolds [20], proposed a Lagrangian stochastic modeling approach for pair dispersion in highly inhomogeneous turbulent flows such as canopy flows.Their approach is consistent with the necessary conditions introduced by Thompson [76], and their model was shown to agree with the experimental results of a scalar released from a known source.Their model also uses the hypothesis that in highly inhomogeneous flows the correlation function for relative velocities is "short".In our measurements of the correlation function (Fig. 8a), we show explicitly at what time lag values this hypothesis is valid.In particular, ρ v || = 0 is reached for all tested initial separations at approximately τ ≈ τ 0 (Fig. 8a).This .Lastly, it is important to regard our interpretation of the results using Richardson's turbulent diffusivity framework.Indeed, turbulent pair dispersion is an intermittent process [10,12,58,65,72], and this characteristic cannot be described using the diffusivity model since it precludes longer correlation times.
Nevertheless, as this is the first empirical investigation of the topic in a canopy flow, it is of particular importance to lay out the leading mechanisms that govern the phenomenon, which include the growth of the kinetic energy with the scale.Since the diffusivity framework lends itself somewhat naturally to the description of this phenomenology, it is the one that was chosen.With this, the treatment of more refined characteristics of pair dispersion at small scales is left for future research.To that should be added the treatment of pair dispersion over larger scales, which could not have been examined using the presently available dataset.

Figure 1 :
Figure 1: Sketches representing the area of measurement.(a) an isometric view of the coordinate system, several roughness elements, and a cutout laser beam illuminating a single sub-volume with tracer particles.(b) a top view of the measurement volume showing the 4 horizontal sub-volume positions.

Figure 2 :
Figure 2: (a) Ratio of the shear length scale and the Kolmogorov length scale for the various sub-volumes plotted as a function of height.(b) The Lagrangian velocity decorrelation timescale, T x , normalized using the mean shear rate plotted as a function of height.Data is shown for the Re ∞ = 26 × 10 3 case.

Figure 3 :
Figure 3: The various components of the pair dispersion tensor are shown for two initial separation values and at the five height groups used.Different shapes correspond to different components of ∆ ij .Lines with the same shape come from each of the four horizontal sub-volume locations and thus represent the horizontal variability of the statistics.

Figure 5 :
Figure 5: Trajectories of pair dispersion anisotropy on the x b y b plane.Two datasets are shown: for five heights with 15ηr 0 < 20η (left) and for four initial separation values with 0.9H < z < 1.1H (right).The beginning of each trajectory, i.e. at time zero, is marked by a black circle, where the trajectories from which the trajectories evolve with time up to 7τ η .Data is for the Re ∞ = 2.6 × 10 4 case, and horizontally averaged across all sub-volumes.

Figure 6 :
Figure 6: Probability distributions of the variance of change in separation distance conditioned on the initial orientation and divided by the unconditioned value.Data is taken at different times and all subvolumes for pairs with 30η < r 0 < 50η.

Figure 7 :
Figure 7: (a) Variance of the change in pair separation normalized with the second moment of separation velocities and the Batchelor timescale.The dashed line corresponds to the Ballistic regime, Eq. (13).(b) Same data but normalized using the time lag, and the dashed line corresponds to Eq. (15).

Figure 8 :
Figure 8: (a) Autocorrelation functions of the rate of separation v || , plotted against time normalized with τ 0 .(b) root mean squared separation distance between trajectories, normalized with the initial separation, and plotted against time lag normalized with the Kolmogorov timescale.

Figure 9 :
Figure 9: The diffusivity of pair dispersion defined as in Eq. (17) for various r 0 cases, plotted against the RMS of the separation distance; (a) data normalized by dissipation scales and (b) by inertial range scaling.Dashed lines correspond to the 4/3s law scaling, Eq. (2).
suggests that the concentration variance can be estimated from Cohen and Reynolds dispersion modelfor times τ > τ v || ≈ τ b =