1. 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 (Reference Batchelor1952), 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 (Raupach & Thom Reference Raupach and Thom1981; Finnigan Reference Finnigan2000; Britter & Hanna Reference Britter and Hanna2003; Belcher Reference Belcher2005; Brunet Reference Brunet2020). In neutrally stable conditions, canopy flows exhibit a distinctive feature characterized by the direct interaction between the 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 (
$\partial U / \partial z$, where U is the mean streamwise velocity component and z is the vertical cordinate), giving rise to coherent Kelvin–Helmholtz structures akin to a mixing-layer analogy (Raupach, Finnigan & Brunet Reference Raupach, Finnigan and Brunet1996). 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, Finnigan & Raupach Reference Brunet, Finnigan and Raupach1994). Interestingly, Lagrangian velocity increments have been observed to follow a Gaussian distribution (Shnapp et al. Reference Shnapp, Bohbot-Raviv, Liberzon and Fattal2020). Additionally, turbulence in canopy flows is typically generated through two interdependent mechanisms: (I) shear production near 
$z \approx H$, and (II) production of smaller-scale turbulence within the wakes of the obstacles within the canopy (
$z< H$) (Finnigan Reference Finnigan2000).
A common framework for modelling dispersion in canopy flows is through Lagrangian stochastic models (DePaul & Sheih Reference DePaul and Sheih1986; Raupach Reference Raupach1987; Baldocchi Reference Baldocchi1997; Katul et al. Reference Katul, Oren, Ellsworth, Hsieh and Phillips1997; Reynolds Reference Reynolds1998; Aylor & Flesch Reference Aylor and Flesch2001; Poggi, Katul & Albertson Reference Poggi, Katul and Albertson2006; Wilson et al. Reference Wilson, Yee, Ek and d'Amours2009; Duman et al. Reference Duman, Trakhtenbrot, Poggi, Cassiani and Katul2016; Bailey, Stoll & Pardyjak Reference Bailey, Stoll and Pardyjak2018; Fattal et al. Reference Fattal, David-Saroussi, Klausner and Buchman2021, Reference Fattal, David-Saroussi, Buchman, Tas and Klausner2023). In such models, Lagrangian fluid particles are advanced through the flow field using random walks that simulate the turbulent flow; this allows estimates of the 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 (Raupach & Thom Reference Raupach and Thom1981; Raupach Reference Raupach1987), 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 (Thomson Reference Thomson1987), the effects of coherent structures (Raupach et al. Reference Raupach, Finnigan and Brunet1996; Ghisalberti & Nepf Reference Ghisalberti and Nepf2002), non-Gaussian velocity PDFs (Pope & Chen Reference Pope and Chen1990), the parallel contributions of wake and shear production (Shnapp et al. Reference Shnapp, Bohbot-Raviv, Liberzon and Fattal2020) or the mechanical diffusion (Nepf Reference Nepf1999). 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 (Batchelor Reference Batchelor1952). 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 wind-tunnel set-up, as we reported in Shnapp et al. (Reference Shnapp, Shapira, Peri, Bohbot-Raviv, Fattal and Liberzon2019). We focused in the past on single-particle statistics of this canopy flow, revealed short decorrelation Lagrangian time scales (Shnapp et al. Reference Shnapp, Bohbot-Raviv, Liberzon and Fattal2020), and characterized the unique features of Lagrangian intermittency in this flow (Shnapp Reference Shnapp2021). 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 modelling 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.
1.1. Background on pair dispersion in isotropic and anisotropic turbulence
Turbulent pair dispersion was introduced by Richardson (Reference Richardson1926) and, since then, it has been studied extensively, e.g. as reviewed in Salazar & Collins (Reference Salazar and Collins2009), Falkovich, Gawdzki & Vergassola (Reference Falkovich, Gawdzki and Vergassola2001) and Sawford (Reference Sawford2001). Pair dispersion describes the statistics of vector 
$\boldsymbol {r}$ connecting two, initially close, Lagrangian fluid particles
\begin{equation} \boldsymbol{r} = \boldsymbol{x}^{(1)} - \boldsymbol{x}^{(2)} ; \quad r \equiv |\boldsymbol{r}| ,\end{equation}
where 
$\boldsymbol {x}^{(i)}$ is the three-dimensional (3-D) 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 (Batchelor Reference Batchelor1952), during which particles’ relative motion is co-linear. (II) At larger time scales, Richardson (Reference Richardson1926) suggested that statistics of 
$r$ can be modelled as a diffusive process with a scale-dependent diffusivity parameter, 
$K$
\begin{equation} K \equiv \frac{1}{2} \frac{\partial \langle {r^2} \rangle}{\partial t} \quad \text{with } K \propto \langle r^2 \rangle^{2/3} \sim r^{4/3} ,\end{equation}
where 
$\langle {\cdot } \rangle$ is an ensemble average of many particle pairs. The scaling of 
$K$, commonly termed ‘the four-thirds law’, leads to a super-diffusive growth of the separation distance with 
$\langle {r^2} \rangle \sim \tau ^{3}$ (Monin & Yaglom Reference Monin and Yaglom1972). (III) At later times, when the two particles have separated farther away than the integral turbulence scale, i.e. for 
$r \gg L$, the distance 
$r$ is expected to be diffusive with a constant, or Taylor, diffusivity (Taylor Reference Taylor1922)
\begin{equation} K = T_L \tilde{u}^2,\end{equation}
where 
$T_L$ is the Lagrangian time scale and 
$\tilde {u}$ is the root-mean-square (r.m.s.) 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 (Salazar & Collins Reference Salazar and Collins2009; Monin & Yaglom Reference Monin and Yaglom1972):
\begin{equation} \langle {(r-r_0)^2} \rangle = \begin{cases} \frac{11}{3}C_2 (\epsilon r_0)^{2/3} \tau^2 & \text{for } \tau \ll \tau_b \\ g \epsilon \tau^3 & \text{for } \tau_b \ll \tau \ll T_L \\ 2 \tilde{u} T_L \tau & \text{for } T_L \ll \tau \end{cases},\end{equation}
where 
$C_2\approx 2.1$ (Sreenivasan Reference Sreenivasan1995) is the so-called Kolmogorov constant, 
$\epsilon$ is the mean rate of the turbulent kinetic energy dissipation, 
$g\approx 0.5$ (Ott & Mann Reference Ott and Mann2000) is the so-called universal Richardson constant and the Batchelor time scale is defined as
\begin{equation} \tau_b = \left( \frac{r_0^2}{\epsilon} \right)^{1/3} ,\end{equation}
with 
$\tau$ the time elapsed since the pair began their separation, starting from 
$r_0$.
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/\eta$ (where 
$\eta =(\nu ^3/\epsilon )^{1/4}$ is the Kolmogorov dissipation length scale). Verification tests of (1.4) in Yeung (Reference Yeung1994), Ott & Mann (Reference Ott and Mann2000), Boffetta & Sokolov (Reference Boffetta and Sokolov2002a), Biferale et al. (Reference Biferale, Boffetta, Celani, Devenish, Lanotte and Toschi2005), Berg et al. (Reference Berg, Lüthi, Mann and Ott2006), Bourgoin et al. (Reference Bourgoin, Ouellette, Xu, Berg and Bodenschatz2006), Ouellette et al. (Reference Ouellette, Xu, Bourgoin and Bodenschatz2006b), Elsinga, Ishihara & Hunt (Reference Elsinga, Ishihara and Hunt2021) and Shnapp et al. (Reference Shnapp, Brizzolara, Neamtu-Halic, Gambino and Holzner2023) 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 
$\tau \gg \tau _b$, and high Reynolds number flows (Berg et al. Reference Berg, Lüthi, Mann and Ott2006; Bourgoin et al. Reference Bourgoin, Ouellette, Xu, Berg and Bodenschatz2006; Ouellette et al. Reference Ouellette, Xu, Bourgoin and Bodenschatz2006a). 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/\eta$ (Biferale et al. Reference Biferale, Boffetta, Celani, Devenish, Lanotte and Toschi2005). Another issue is the strong intermittency of pair dispersion (Bitane, Homann & Bec Reference Bitane, Homann and Bec2012; Scatamacchia, Biferale & Toschi Reference Scatamacchia, Biferale and Toschi2012; Biferale et al. Reference Biferale, Lanotte, Scatamacchia and Toschi2014; Thalabard, Krstulovic & Bec Reference Thalabard, Krstulovic and Bec2014; Shnapp & Liberzon Reference Shnapp and Liberzon2018), presumably due to Lagrangian statistics associated with long-time correlations (Boffetta & Sokolov Reference Boffetta and Sokolov2002b). In particular, recent observations have shown that the scaling of 
$\langle {r^2} \rangle$ depends on the initial conditions, namely on the initial separation (Elsinga et al. Reference Elsinga, Ishihara and Hunt2021; Tan & Ni Reference Tan and Ni2022; Shnapp et al. Reference Shnapp, Brizzolara, Neamtu-Halic, Gambino and Holzner2023) and the initial separation velocity (Shnapp & Liberzon Reference Shnapp and Liberzon2018), while it is also known that particles may retain their separation distance for very long times (Scatamacchia et al. Reference Scatamacchia, Biferale and Toschi2012). Thus, it was suggested that Richardson's diffusive approach relies on an assumption of short Lagrangian correlation times (Bitane et al. Reference Bitane, Homann and Bec2012), 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 (Sokolov, Klafter & Blumen Reference Sokolov, Klafter and Blumen2000; Rast & Pinton Reference Rast and Pinton2011; Thalabard et al. Reference Thalabard, Krstulovic and Bec2014; Bourgoin Reference Bourgoin2015). Bitane et al. (Reference Bitane, Homann and Bec2012) proposed an alternative to 
$\tau _b$ which renders the process self-similar with respect to 
$r_0$; this framework is discussed in more detail below.
While (1.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 
$\boldsymbol {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 (Batchelor Reference Batchelor1952)
\begin{equation} \varDelta_{ij}(\tau) = \langle { (r_i - r_{i,0} ) {\cdot} (r_j - r_{j,0})} \rangle ;\end{equation}
notably, 
$\mathrm {tr}(\varDelta _{ij}) = \langle {(r-r_0)^2} \rangle$, and in HIT 
$\varDelta _{ij}$ is diagonal. Furthermore, pair dispersion in anisotropic flows can depend on the initial orientation of 
$\boldsymbol {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 & Yeung (Reference Shen and Yeung1997) established that the presence of a mean shear has a strong effect on pair dispersion: it makes the evolution of 
$\varDelta _{ij}$ anisotropic and with non-zero non-diagonal components, and it can render apparent super-diffusivity at long times (
$\varDelta _{xx}\sim \tau ^3$) irrespective of the argument of scale dependence. More recently, Pitton et al. (Reference Pitton, Marchioli, Lavezzo, Soldati and Toschi2012) and Polanco et al. (Reference Polanco, Vinkovic, Stelzenmuller, Mordant and Bourgoin2018) 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 direct numerical simulation of a square duct flow, Sharma & Phares (Reference Sharma and Phares2006) reported that particles remained trapped in secondary flow regions which affected the power law behaviour of 
$r(\tau )$. Celani et al. (Reference Celani, Cencini, Vergassola, Villermaux and Vincenzi2005) 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 fluctuation dominance or mean shear dominance. A cross-over between these two regimes was predicted to occur at a critical time scale that is proportional to the mean shear 
$\tau _c\sim ({{\rm d}U}/{{\rm d}z})^{-1}$. Recent studies further examined buoyancy-driven flows (Schumacher Reference Schumacher2008; Ni & Xia Reference Ni and Xia2013; Liot et al. Reference Liot, Martin-Calle, Gay, Salort, Chillà and Bourgoin2019) and pair dispersion of inertial particles (Pitton et al. Reference Pitton, Marchioli, Lavezzo, Soldati and Toschi2012) or in other complex flows, such as in the stroke of a swimming jellyfish (Kim & Chamorro Reference Kim and Chamorro2019).
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.
2. Methods
2.1. A 3-D particle tracking wind-tunnel experiment
We study pair dispersion in a canopy flow using the results of a 3-D particle tracking velocimetry experiment (3D-PTV Malik, Dracos & Papantoniou Reference Malik, Dracos and Papantoniou1993; Mass, Gruen & Papantoniou Reference Mass, Gruen and Papantoniou1993; Dracos Reference Dracos1996). In the experiment, we obtained flow tracer trajectories in a wind-tunnel canopy flow model using an extended, real-time image processing system (Shnapp et al. Reference Shnapp, Shapira, Peri, Bohbot-Raviv, Fattal and Liberzon2019). The wind-tunnel laboratory, situated at the Israel Institute for Biological Research (Ness Ziona, Israel), features a 14 metre long, open circuit suction wind tunnel with a 
$2\times 2$ squared metre cross-sectional area, that is fit for conducting experiments that mimic turbulent flows in the atmospheric surface layer (Bohbot-Raviv et al. Reference Bohbot-Raviv, Shnapp, Liberzon, Babin, Hotoveli, Shick and Fattal2017). The canopy flow was modelled by placing flat rectangular plates along the entire bottom wall of the test section (figure 1d). Our study used an inhomogeneous canopy layer, constructed with two types of flat plates with heights either 
$H$ or 
$\frac {1}{2}H$, and width 
$\frac {1}{2}H$, where 
$H=100$ mm; the two types of plates were positioned in consecutive rows and at a staggered construction, see figure 1. The canopy frontal area density, defined as 
$\varLambda _f = A_f /A_T$, (
$A_f$ being the element frontal area, and 
$A_T$ the lot area of the canopy) is 
$\frac {1}{2}$, which positions our canopy between the ‘dense’ and ‘sparse’ categories (Brunet Reference Brunet2020). Data were gathered at two levels of the free-stream velocity, 
$U_\infty = 2.5$ and 
$4 \ {\rm m} \ {\rm s}^{-1}$, being the free-stream mean velocity measured with a sonic anemometer at the centre of the wind-tunnel cross-section. These values correspond to Reynolds numbers 
$Re_\infty \equiv {U_\infty H}/{\nu } \approx 16\times 10^3$ and 
$26\times 10^3$, where 
$\nu$ is the kinematic viscosity of the air.
Figure 1. The canopy layer model used in the experiment. (a) An isometric view schematic sketch of the coordinate system, several roughness elements and a cutout laser beam illuminating a single sub-volume with tracer particles. (b) A top view schematic sketch of the measurement volume showing the 4 horizontal sub-volume positions. (c) An image of the canopy roughness elements with the laser beam and seeding particles seen passing through it. (d) A part of the wind-tunnel test section with the roughness elements placed on the bottom wall; the four camera system is seen outside of the test section, directed at the measurement region.
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 Shnapp et al. (Reference Shnapp, Shapira, Peri, Bohbot-Raviv, Fattal and Liberzon2019), so only brief information shall be repeated here for completeness. The flow was seeded with hollow glass spheres with an estimated Stokes number of 
$\overline {St}\approx 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 particles were added to the air flow inside the test section, far upstream of the measurement region through four nuzzles connected to a fluidized bed reactor. The particles were then carried by the air flow in the wind tunnel down to the measurement region, where they were illuminated with a 10 W continuous wave laser beam, as shown in figure 1(a) and 1(c). A set of four 4-megapixel cameras captured images of the tracer particles at a rate of either 500 or 1000 Hz, depending on the height at which data were gathered. The images from each camera were fed into a dedicated FPGA card for online image analysis and particle segmentation on the hardware, and the particle two-dimensional (2-D) positions were saved for post-processing in a compact binary format. This recording system, internally called ‘Blob Recorded’, allows unprecedentedly long 3D-PTV experimental campaigns to be performed due to the achievable 
$O(10^4)$ data compression ratio, see Shnapp et al. (Reference Shnapp, Shapira, Peri, Bohbot-Raviv, Fattal and Liberzon2019). Following that, we used the open source software OpenPTV (OpenPTV Consortium 2004) to perform camera calibration, stereo matching and tracking. Finally, particle positions were smoothed and the velocities and accelerations were calculated by local spline fitting following the method introduce in Luthi, Tsinober & Kinzelbach (Reference Luthi, Tsinober and Kinzelbach2005). To ensure good trackability of the particles in the high-speed wind-tunnel air flow we used rather low seeding values of the order of 20 particles per image; together with very long recording durations, approximately 15 minutes long in each location, we ensured that the statistics we gathered represent the flow well. The imaging resolution we used was 50 
$\mathrm {\mu }$m per pixel, and the region of measurement was divided into sub-volumes, approximately 
$8\times 5\times 4\ {\rm mm}^3$ in size (full description is given in § 2.3). The data analysis was performed using the open source Flowtracks package (Meller & Liberzon Reference Meller and Liberzon2016). 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.
2.2. The canopy flow model
The canopy model in the experiment had a flow structure with characteristics that resemble those typically observed in neutrally buoyant canopies (Finnigan Reference Finnigan2000), however, with certain unique features that were the result of the double-height construction we used. The double-averaged mean velocity profile demonstrates an inflection point slightly above the roughness elements. Our previous work with this canopy construction, based on particle image velocimetry measurements, showed a smaller shear length at the top of the canopy as compared with homogeneous canopy layers at similar densities and that, unlike homogeneous canopies, it had a second inflection point right above the lower roughness elements (Shig et al. Reference Shig, Babin, Shnapp, Fattal, Liberzon and Bohbot-Raviv2023a). The velocity component distributions were anisotropic, they become broader with 
$z$, and a pronounced positive skewness was present in the distributions of the streamwise velocity component (Shnapp et al. Reference Shnapp, Shapira, Peri, Bohbot-Raviv, Fattal and Liberzon2019); this skewness is a hallmark of the sweep events that characterize the roughness layer (
$z<1.5H$) in canopy flows (Brunet et al. Reference Brunet, Finnigan and Raupach1994; Raupach et al. Reference Raupach, Finnigan and Brunet1996). In addition to that, a second peak of increased sweep contribution above the lower set of elements appeared due to the double-height canopy construction (Shig et al. Reference Shig, Babin, Shnapp, Fattal, Liberzon and Bohbot-Raviv2023b).
From the Lagrangian statistics point of view, our previous analysis of 3D-PTV measurements showed that the canopy is characterized by an unexpectedly short Lagrangian de-correlation time scale (Shnapp et al. Reference Shnapp, Bohbot-Raviv, Liberzon and Fattal2020). This feature leads to PDFs of velocity increments that are of a Gaussian shape despite the marked non-Gaussianity of the absolute velocity components. In addition to that, single-particle velocity statistics were found to exhibit both the characteristic small-scale Lagrangian intermittency (Arnèodo et al. Reference Arnèodo2008), and a so-called large-scale intermittency due to the sweep–ejection cycle of the canopy flow (Shnapp Reference Shnapp2021).
2.3. 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 figure 1(a). 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 labelled alphabetically 
$a$, 
$b$, 
$c$ and 
$d$, as shown in figure 1(b). At each position, 
$a\unicode{x2013}d$, we used 5 vertical slabs of thickness 
$0.2 H$ to define the sub-volumes; the vertical slab positions are labelled numerically 1–5 and correspond to heights 
$z/H \in \{0.5\unicode{x2013}0.7, 0.7\unicode{x2013}0.9, 0.9\unicode{x2013}1.1, 1.1\unicode{x2013}1.3, 1.3\unicode{x2013}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
\begin{equation} \langle {A} \rangle \equiv \frac{1}{N}\sum_{i=0}^N A_i(\boldsymbol{x}_{\mathbb{V}},t),\end{equation}
where 
$\mathbb {V}$ is a tag pertaining to different particle ensembles. In our recent work (Shnapp et al. Reference Shnapp, Bohbot-Raviv, Liberzon and Fattal2020), we explored single-particle statistics using the same dataset. There, the sub-volume-averaged quantity 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 the random turbulent forcing outweighed the combined effects of flow inhomogeneity terms in the stochastic particle equation of motion. Specifically, this was achieved through the condition
\begin{equation} \left\langle { \frac{1}{2} \frac{|C_0 \epsilon R_{ij}^{{-}1} u_j'|}{|\phi_i /g|} } \right\rangle \ll 1,\end{equation}
where 
$C_0$ is the so-called Lagrangian structure function coefficient, 
$R_{ij}=\langle {u'_i \, u_j'} \rangle$ is the Reynolds stress tensor, 
$u'_j$ is the r.m.s. of the 
$j$ velocity component and 
$\phi _i/g$ is a vector representing a sum of drift terms of the Lagrangian velocity PDF that allows us to account for flow inhomogeneity (Thomson Reference Thomson1987). In the present work, we capitalize on the condition (2.2), 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 Shnapp et al. (Reference Shnapp, Bohbot-Raviv, Liberzon and Fattal2020), turbulent parameters of our flow were calculated where they were discussed and analysed in depth. In particular, we calculated 
$\epsilon$, the turbulence dissipation length and time scales 
$\eta$, 
$\tau _\eta$, the Lagrangian decorrelation time scale 
$T_{i}$ and the Taylor Reynolds number 
$Re_\lambda$. The values of these parameters are tabulated for each sub-volume in the Appendix, and they shall be used in the analysis below.
3. Local isotropy of pair dispersion
3.1. Scaling analysis
As discussed in § 1, numerical simulations (Shen & Yeung Reference Shen and Yeung1997; Sharma & Phares Reference Sharma and Phares2006; Pitton et al. Reference Pitton, Marchioli, Lavezzo, Soldati and Toschi2012; Polanco et al. Reference Polanco, Vinkovic, Stelzenmuller, Mordant and Bourgoin2018; Polanco Reference Polanco2019) confirmed theoretical predictions (Monin & Yaglom Reference Monin and Yaglom1972) that the mean shear, inhomogeneity and/or anisotropy in the flow cause anisotropic pair dispersion. This 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$. Empirical evidence shows that turbulent velocity fluctuations are characterized by a tendency to recover isotropy at small scales, a phenomenon dubbed the return to isotropy (Lumley & Newman Reference Lumley and Newman1977). The return to isotropy at smaller scales was also observed in canopy flows (Stiperski, Katul & Calaf Reference Stiperski, Katul and Calaf2021). 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 Celani et al. (Reference Celani, Cencini, Vergassola, Villermaux and Vincenzi2005). 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 (Finnigan Reference Finnigan2000): (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 with form drag due to the low velocity near the ground for sufficiently dense canopies. Because shear production occurs due to the interaction of the double-averaged mean shear and the Reynolds stress, 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 elements (Finnigan Reference Finnigan2000). Since the local orientation of most straining directions changes in space, we expect that the anisotropy in the 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
\begin{equation} L_\varGamma = \left( \frac{\epsilon}{\varGamma^3} \right)^{1/2},\end{equation}
where 
$\varGamma = {{\rm d}U}/{{\rm d}z}$ is the mean velocity gradient. Indeed, analysing the scale-by-scale turbulent kinetic energy budget in turbulent shear flows, Casciola et al. (Reference Casciola, Gualtieri, Benzi and Piva2003) showed that the production by shear continues down to 
$L_\varGamma$, whereas for scales 
$r< L_\varGamma$ the turbulent kinetic energy transport becomes more dominant. For our analysis, we parameterize the shear length using a global, canopy-wide, measure of the mean velocity gradient 
${-\varGamma _z = {U(z=1.5H)}/{1.5H} \approx 4.75\ {\rm s}^{-1}}$. Figure 2(a) shows that 
$L_\varGamma /\eta$ increased from roughly 100 inside the canopy and up to roughly 150 above it.
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 time scale, 
$T_x$, normalized using the mean shear rate plotted as a function of height. Data are shown for the 
$Re_\infty =26\times 10^3$ case.
Another critical dimension of pair dispersion is time. Celani et al. (Reference Celani, Cencini, Vergassola, Villermaux and Vincenzi2005) 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 time scale for this transition is proportional to the inverse mean shear rate, 
$\tau _c \sim \varGamma ^{-1}$. Figure 2(b) shows the Lagrangian velocity decorrelation time scale for the streamwise velocity component normalized by 
$\varGamma$. Throughout our measurement region, we observe that 
$T_x {\cdot } \varGamma <1$.
Due to the finite measurement region, our work focuses on time scales up to 
$\tau \approx T_x$. In addition, the typical separations between particles that we consider in our work reach up to roughly 
$r=100\eta$, namely, mostly values with 
$r< L_\varGamma$. Thus, the dimensions of our measurement region and the time scales on which we focus are smaller than the scales imposed by the mean shear. For this reason, we expect to see only weak anisotropy in pair-dispersion statistics due to the local isotropy of the flow.
3.2. Observation of weak anisotropy
The various components of the pair-dispersion tensor, 
$\varDelta _{ij}$, across the whole region of measurements are shown in figure 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 separations in different directions. This behaviour is unlike what is expected in a mean shear-driven separation scenario as, for example, in homogeneously sheared turbulence the separation velocity in the streamwise direction grows as the separation along the mean velocity gradient direction increases. Thus, the behaviour observed for 
$\varDelta _{ij}$ is consistent with the notion presented in § 3.1 that the mean shear is not expected to significantly affect pair dispersion in our measurements.
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 
$\varDelta _{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.
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:
\begin{equation} I_{ij} \equiv \frac{\varDelta_{ij}}{\mathrm{tr}(\varDelta_{ij})} -\frac{1}{3}\delta_{ij},\end{equation}
where 
$\delta _{ij}$ is the identity matrix. The tensor 
$I_{ij}$ shows the statistics of the components of 
$\boldsymbol {{\rm \Delta} r}$ with respect to their norms, essentially highlighting anisotropy. Let us note that 
$\mathrm {tr}(I_{ij})=\sum _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 
$\frac {1}{3}\langle {(r-r_0)^2} \rangle$), while for 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}=\frac {2}{3}$ while 
$I_{yy}=I_{zz}=-\frac {1}{3}$.
Histograms of the diagonal components, 
$I_{ii}(\tau )$ are shown in figure 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 positions, both 
$Re_\infty$, and across the time range 
$\tau <6\tau _\eta$, during which 
$I_{ii}$ did not change appreciably.
Figure 4. Diagonal terms of the pair-dispersion tensor normalized by its trace minus one third. Data are shown for sub-volumes at three heights and for two levels of 
$r_0$.
The disparity between the separation components decreases with height, being the largest inside the canopy and smallest above it. Furthermore, the increase in separation was the slowest in the vertical direction and the fastest in the spanwise direction, whereas the increase of the streamwise separation was approximately at isotropic values 
$(I_{xx} \approx 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 channelling of the flow in between the roughness obstacles, which were in a staggered configuration, and not necessarily due to the effects of the 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 (Shen & Yeung Reference Shen and Yeung1997; Pitton et al. Reference Pitton, Marchioli, Lavezzo, Soldati and Toschi2012; Polanco et al. Reference Polanco, Vinkovic, Stelzenmuller, Mordant and Bourgoin2018).
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 (Lumley & Newman Reference Lumley and Newman1977).
The increase 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. (Reference Stiperski, Katul and Calaf2021), which used a projection of the invariants of the normalized Reynolds stress tensor on a 2-D 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
\begin{equation} \begin{array}{c@{}}
x_b = \lambda_1 - \lambda_2 + \dfrac{1}{2}(3\lambda_3 +
1)\\ y_b = \dfrac{\sqrt{3}}{2} (3\lambda_3 + 1)
\end{array},\end{equation}
where 
$\lambda _1$, 
$\lambda _2$ and 
$\lambda _3$ are the smallest, intermediate and largest eigenvalues of 
$I_{ij}$, respectively. Equation (3.3) maps the eigenvalues to a planar equilateral triangle. As seen in figure 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 us 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 
$\varGamma ^{-1}$.
Figure 5. Trajectories of pair-dispersion anisotropy on the 
$x_b$–
$y_b$ plane. Two datasets are shown: for five heights with 
$15\eta r_0<20\eta$ (a) and for four initial separation values with 
$0.9H< z<1.1H$ (b). The beginning of each trajectory, i.e. at time zero, is marked by a black circle, from which the trajectories evolve with time up to 7
$\tau _\eta$. Data are for the 
${Re}_\infty =2.6\times 10^4$ case, and horizontally averaged across all sub-volumes.
3.3. Bias due to initial orientation
As shown in Shen & Yeung (Reference Shen and Yeung1997), Pitton et al. (Reference Pitton, Marchioli, Lavezzo, Soldati and Toschi2012) and Polanco et al. (Reference Polanco, Vinkovic, Stelzenmuller, Mordant and Bourgoin2018), 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 
$\boldsymbol {u}= (\varGamma z, 0, 0)$ and considering particles with an initial separation 
$\boldsymbol {r}_0=(r_{0,x}, r_{0,y}, r_{0,z})$, then the separation vector is 
$\boldsymbol {r}(\tau ) = (r_{0,x} + \tau \varGamma 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 
$\boldsymbol {r}_o$ is aligned with the mean shear direction, the bias is by a factor of 
$\varGamma r_{0}$.
The relative importance of this bias in a turbulent flow can be estimated by comparing 
$\varGamma r_0$ with the spatial fluctuating velocity increments 
$\delta _{r_0} u$. If we apply the Kolmogorov dimensional scaling (Kolmogorov Reference Kolmogorov1941) to parameterize 
$\delta _{r_0} u \sim (\epsilon r_0)^{1/3}$, we can construct a dimensionless group 
$S = {\varGamma r_0}/{\delta _{r_0} u} = {\varGamma r_0^{2/3}}/{\epsilon ^{1/3}} = ({r_0}/{L_\varGamma })^{2/3}$ that quantifies the importance of the shear bias for particles with 
$r$ in the inertial range. When 
$S \gg 1$ we expect that shear will bias pair separation based on their initial orientation whereas when 
$S \ll 1$ it will not. Also, if an inertial range exists, this bias grows stronger as 
$r_0^{2/3}$ for smaller initial separations, which is in qualitative agreement with the observations of Shen & Yeung (Reference Shen and Yeung1997) and Polanco et al. (Reference Polanco, Vinkovic, Stelzenmuller, Mordant and Bourgoin2018). Also, in the dissipation range, the velocity field is presumably smooth and the typical size of velocity gradients is 
${1}/{\tau _\eta } = ({\epsilon }/{\nu })^{1/2}$ (where 
$\tau _\eta$ is the dissipation time scale). Therefore, the dissipation scaling for velocity increments in the dissipation range is 
$\delta _{r_0} u \sim r_0 ({\epsilon }/{\nu })^{1/2}$ so that the dimensionless parameter becomes 
$S = \varGamma ({\epsilon }/{\nu })^{-1/2}$. To summarize, we define
\begin{equation} S = \begin{cases} \varGamma \left( \dfrac{\epsilon}{\nu} \right)^{{-}1/2} & \text{if } r_0\ll \eta\\ \left( \dfrac{r_0}{L_\varGamma} \right)^{2/3} & \text{if } r_0\gg \eta \end{cases}, \end{equation}
and when 
$S\gg 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 
$\langle {\delta r^2} \rangle$ conditioned on the initial orientation of 
$\boldsymbol {r}_0$. We divided the pairs into sub-samples based on the condition that the angle between 
$\boldsymbol {r}_0$ and each of the coordinate axes was lower than 
$25^\circ$. In figure 6, probability distributions are shown for the conditional variance normalized by the variance of all the pairs. In this case, data were taken at all available times, all sub-volumes, and with 
$30\eta < r_0 < 50 \eta$. The distributions show that pairs with initial separation aligned vertically (i.e. with 
$\hat {z}$) typically separated faster than the average, whereas pairs with initial separation aligned with the spanwise (i.e. 
$\hat {y}$) typically separated slower than the average. Considering the previous observation that separations are fastest in the spanwise direction (figure 4), the bias observed here for vertically oriented pairs may suggest that the channelling effect varies with height. This is also consistent with the observation of the strongest anisotropy inside the canopy layer.
Figure 6. Probability distributions of the variance of change in separation distance conditioned on the initial orientation and divided by the unconditioned value. Data are taken at different times and all sub-volumes for pairs with 
$30\eta < r_0 < 50\eta$.
From figure 2 it is seen that for the present case, we have a maximum value of 
$S\approx 0.38$, which suggests that the bias should be quite weak. This is validated in figure 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 (Pitton et al. Reference Pitton, Marchioli, Lavezzo, Soldati and Toschi2012; Polanco et al. Reference Polanco, Vinkovic, Stelzenmuller, Mordant and Bourgoin2018) 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. (Reference Celani, Cencini, Vergassola, Villermaux and Vincenzi2005). These results reveal the important role of 
$L_\varGamma$ and 
$\varGamma$ 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.
4. Regimes of pair dispersion and a scale-dependent diffusivity
Following the observation in § 3 that pair dispersion is only weakly anisotropic in our measurements, we focus next on statistics of the change of the separation distance, 
$r=|\boldsymbol {r}|$, with time. In particular, we analyse the scaling of the change in separation distance 
$\langle {(r-r_0)^2} \rangle \propto \tau ^\beta$ (Biferale et al. Reference Biferale, Boffetta, Celani, Devenish, Lanotte and Toschi2005; Berg et al. Reference Berg, Lüthi, Mann and Ott2006; Ouellette et al. Reference Ouellette, Xu, Bourgoin and Bodenschatz2006b; Salazar & Collins Reference Salazar and Collins2009; Shnapp & Liberzon Reference Shnapp and Liberzon2018). It is also noted that, formally, due to the inhomogeneity of our canopy flow, 
$\beta$ is essentially an ‘effective scaling’. Therefore, the goal of this section is to compare the empirical data with the different regimes of (1.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 Biferale et al. (Reference Biferale, Bodenschatz, Cencini, Lanotte, Ouellette, Toschi and Xu2008). 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 
$\beta$. 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 \tau _\eta$. Statistics were calculated only for 
$\tau \leq 5 \tau _\eta$ (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.
4.1. The ballistic regime
We first verify that pair dispersion at short times follows a ballistic propagation, as was suggested by Batchelor (Reference Batchelor1952). This regime is characterized by a quadratic growth of 
$\langle {{\rm \Delta} r^2} \rangle \propto \tau ^2 (\tau \ll \tau _b)$, corresponding to the first case of (1.4). Figure 7(a) presents pair dispersion from the canopy flow in normalized form according to (1.4)
\begin{equation} \frac{\langle {{\rm \Delta} r^2} \rangle}{S_{LL}(r_0) \tau_b} = \left( \frac{\tau}{\tau_b} \right)^2,\end{equation}
where 
$S_{LL}(r_0) = \langle {(\delta _{r_0} \boldsymbol {u} \boldsymbol {\cdot } {\boldsymbol {r_0}}/{r_0})^2} \rangle$ is the directly measured Eulerian structure function, and the Batchelor time scale 
$\tau _b$. Initially, for 
$\tau < 0.1 \tau _b$, all the curves with different initial separations collapse on a quadratic dashed line, as (1.4) suggests. Similar results were previously obtained in homogeneous flows by Ouellette et al. (Reference Ouellette, Xu, Bourgoin and Bodenschatz2006b), Bitane et al. (Reference Bitane, Homann and Bec2012) and Shnapp & Liberzon (Reference Shnapp and Liberzon2018), and inhomogeneous flows by Polanco et al. (Reference Polanco, Vinkovic, Stelzenmuller, Mordant and Bourgoin2018), Liot et al. (Reference Liot, Martin-Calle, Gay, Salort, Chillà and Bourgoin2019) and Ni & Xia (Reference Ni and Xia2013). 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$.
4.2. 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. (Reference Bitane, Homann and Bec2012) for the case of HIT flow. We denote the separation velocity 
$v_{\|} = {\partial r}/{\partial t}$ and acceleration 
$a_{\|} = {\partial v_{\|}}/{\partial t}$. Then, the separation distance is Taylor expanded to the third order, squared and ensemble averaged, which gives
\begin{equation} \langle {{\rm \Delta} r^2} \rangle = \langle {v_{\|,0}^2} \rangle \tau^2 + \tfrac{1}{2} \langle {a_{\|,0} \, v_{\|,0}} \rangle \tau^3 + O(\tau^4),\end{equation}
where subscript 
$0$ denotes initial values at 
$\tau =0$. This expression, which extends the Ballistic regime, is expected to hold for short times, and, in particular, according to Bitane et al. (Reference Bitane, Homann and Bec2012), it gives the appropriate time scale for the end of the ballistic regime at 
$\tau _0 = -{\langle {v_{\|,0}^2} \rangle }/{\langle {a_{\|,0} \, v_{\|,0}} \rangle }$. Note that the minus sign is used since in 3-D turbulence 
$\langle {a_{\|,0} \, v_{\|,0}} \rangle < 0$ (Ott & Mann Reference Ott and Mann2000), as was confirmed for our data (not shown for brevity). Therefore, the above Taylor expansion can be rewritten in dimensionless form as
\begin{equation} \frac{\langle {{\rm \Delta} r^2} \rangle}{\langle v_{\|,0}^2 \rangle \tau^2} = 1 - \frac{1}{2} \left( \frac{\tau}{\tau_0} \right) + O(\tau^2).\end{equation}
In figure 7(b) we show the empirical data of 
$\langle {{\rm \Delta} r^2} \rangle$ for the different 
$r_0$ cases, normalized according to the left-hand side of (4.3), and compare with the right-hand side shown as a dashed black line. For short times, 
$\tau \ll \tau _0$, the figure shows rather good agreement between the empirical data and (4.3), especially in the downward trend of the lines, which means that the above arguments of Bitane et al. (Reference Bitane, Homann and Bec2012) agree with our data. Furthermore, the negative term in (4.3) is responsible for the deviations from a purely quadratic growth as time progresses. Therefore, defining the end of the ballistic regime as a deviation of 
$\langle {{\rm \Delta} r^2} \rangle$ by, say, 10 % from 
$\langle {v_{\|,0}} \rangle \tau ^2$, we obtain that the ballistic regime ends at 
$\tau =0.2 \tau _0$. Indeed, our data from different 
$r_0$ agree with (4.3) (dashed line) roughly up to 
$\tau \approx 0.2 \tau _0$. At longer times, 
$\tau \gtrsim 0.2 \tau _0$, the normalized pair dispersion curves in figure 7(b) begin deviating upwards. The increase in of the curves slope signals an increase in the time scaling of the separation, which can be attributed to the end of the ballistic regime and a transition towards the inertial regime (transition from the first to the second case in (1.4)).
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
\begin{equation} \rho_{v_{\| }}(\tau) = \frac{\langle {v_{\| }(\tau) v_{\| }(0)} \rangle}{ \sqrt{\langle { v_{\| }^2(0)} \rangle \langle {v_{\| }^2(\tau)} \rangle}},\end{equation}
which is shown in figure 8. The autocorrelation 
$\rho _{v_{\|}}$ is plotted against 
${\tau }/{\tau _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 
$\tau \sim \tau _0$ (with some weak dependence on 
$r_0$). This confirms that, at 
$\tau \gtrsim \tau _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 regime.
Figure 8. Autocorrelation functions of the rate of separation 
$v_{\|}$, shown against time normalized with 
$\tau _0$.
4.3. Scale-dependent growth rates and approach to Richardson's 4/3 law
Richardson's classical theory (Richardson Reference Richardson1926) 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 figure 9 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\eta$, 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 (Ott & Mann Reference Ott and Mann2000; Biferale et al. Reference Biferale, Boffetta, Celani, Devenish, Lanotte and Toschi2005; Liot et al. Reference Liot, Martin-Calle, Gay, Salort, Chillà and Bourgoin2019; Elsinga et al. Reference Elsinga, Ishihara and Hunt2021; Tan & Ni Reference Tan and Ni2022), because the separation time scale, 
$\tau _0 \sim r_0^{2/3}$, 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 figure 7(a), an increase of the local slope for the 
$r_0=10\eta$ group can be seen at longer times (
$\tau \ge \tau _0$), indicating an increase in local scaling exponent. In these cases, the pair dispersion is super-diffusive, i.e. 
$\langle {{\rm \Delta} r^2} \rangle \propto \tau ^\beta$ with 
$\beta >2$.
Figure 9. Root mean squared separation distance between trajectories, normalized with the initial separation and plotted against time lag normalized with the Kolmogorov time scale.
With the observations obtained thus far, let us lay out the reasoning to explain the behaviour of pair dispersion in figure 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. figure 8), so that a so-called inertial regime was reached. Notably, even for the largest initial separation considered, 
$r_0=50\eta$, we have still have 
$\tau _0 /T_L\approx 0.96$, so the pairs did not enter yet into the asymptotic diffusive regime (third case in (1.4)). Second, the increase of 
$r$ relative to its initial value depended on 
$r_0$, as shown in figure 9. 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 
$\beta >2$. In contrast, for the largest 
$r_0$ case there was 
${\langle {r^2} \rangle ^{1/2}}/{r_0}\approx 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 
$\beta \approx 1$, as expected in a diffusive process with constant diffusivity.
We support the above picture by estimating diffusivity directly. Following Batchelor (Reference Batchelor1952), we define
\begin{equation} K_r(r) \equiv \frac{1}{2} \frac{\partial \langle {r^2} \rangle}{\partial t}.\end{equation}
In figure 10(a) the diffusivity is plotted against the r.m.s. 
$\tilde {r} \equiv \langle {r^2} \rangle ^{1/2}$, where the axes are normalized using the dissipation scales, 
$\eta$ and 
$\tau _\eta$. 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 Pumir, Shraiman & Chertkov (Reference Pumir, Shraiman and Chertkov2001), is a consequence of the negative skewness of the Eulerian velocity differences, aka the 4/5 law (Monin & Yaglom Reference Monin and Yaglom1972). With increasing time, 
$K_r$ grows. At small times, 
$K_r$ grows in line with the ballistic regime, the end of which, at 
$\tau =0.2\tau _0$, is indicated as circles on the figure. Later on, at 
$\tau \approx \tau _0$, the separation velocity became decorrelated, which is followed by the inertial regime. The time 
$\tau =\tau _0$ is indicated by stars on the figure. The range of our measurement extended to this regime (
$\tau \gtrsim \tau _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 
$\tilde {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 Reference Batchelor1952), according to (1.2) predicts that 
$K_r \propto \tilde {r}^{4/3}$. The prediction (1.2) is plotted in the figure as a dashed black line. The data for the two smallest 
$r_0$ cases lean towards the dashed line for 
$\tau >\tau _0$. This observation agrees with previous results that showed a super-diffusive, 
$t^3$ scaling range for small initial separations in HIT (Elsinga et al. Reference Elsinga, Ishihara and Hunt2021; Tan & Ni Reference Tan and Ni2022).
The Kolmogorov scaling for the diffusivity in the inertial range is 
$K_r \sim (\tilde {r}^{4} \epsilon )^{1/3}$. Therefore, we annotate 
$K_{r_0} \equiv (r_0^{4} \, \epsilon )^{1/3}$. Accordingly, in figure 10(b) the same data are plotted where the axes are normalized using 
$r_0$ and 
$K_{r_0}$. The data from all 
$r_0$ cases are seen to collapse on a single curve under this normalization. The collapse of the data from all 
$r_0$ in figure 10(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 join the 
$K_r = K_{r_0} (\tilde {r}^{4} \epsilon )^{1/3}$ line, which corresponds to Richardson's 4/3 law, (1.2). These observations extend the fundamental phenomenology of turbulent pair dispersion to small-scale separation in highly turbulent anisotropic flows.
5. 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 (Shnapp et al. Reference Shnapp, Shapira, Peri, Bohbot-Raviv, Fattal and Liberzon2019).
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 the 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_\varGamma$ – the scale that determines the cross-over between small-scale local isotropy and large-scale anisotropy (Celani et al. Reference Celani, Cencini, Vergassola, Villermaux and Vincenzi2005). In addition, the time range of our measurement was smaller than 
$\varGamma ^{-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 channelling in between the roughness obstacles. Investigating the pair-dispersion tensor in a framework analogous to the return to isotropy (Lumley & Newman Reference Lumley and Newman1977; Stiperski et al. Reference Stiperski, Katul and Calaf2021) (figure 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 (Stiperski et al. Reference Stiperski, Katul and Calaf2021). This similarity might hint that theory developed in the framework of ‘return to isotropy’ could provide helpful information on pair-dispersion modelling in anisotropic flows; however, this is left for future work.
An important feature observed in our experiments was the smallness of turbulence scales compared with 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_\varGamma / \eta$ and 
$T_L \varGamma$ 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 two initial regimes of pair dispersion, the ballistic and the inertial, 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 (Bourgoin et al. Reference Bourgoin, Ouellette, Xu, Berg and Bodenschatz2006). 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\eta$), in agreement with previous observations in isotropic turbulence (Elsinga et al. Reference Elsinga, Ishihara and Hunt2021; Tan & Ni Reference Tan and Ni2022). 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.
The local isotropy in pair dispersion observed in our canopy flow holds due to the focus of our measurements on relatively small scales, 
$r \ll H$ and 
$r< L_\varGamma$. At longer times we expect the local isotropy to break and that other mechanisms will become dominant. In particular, although inaccessible in our current dataset, we expect three phenomena to be important at longer times and larger separations. First, the strong shear above the canopy (
$z>H$) is expected to make the streamwise component of the separation in this region faster than the vertical and lateral components, thus making pair dispersion at longer times highly anisotropic. This shear dominated regime, characterized by a 
$r\sim t^{3/2}$ scaling (Monin & Yaglom Reference Monin and Yaglom1972), was observed in numerical simulations of homogeneous sheared turbulence and channel flows (Shen & Yeung Reference Shen and Yeung1997; Pitton et al. Reference Pitton, Marchioli, Lavezzo, Soldati and Toschi2012; Polanco et al. Reference Polanco, Vinkovic, Stelzenmuller, Mordant and Bourgoin2018), and is expected to occur also above canopy flows for 
$r\gg L_\varGamma$, according to the theory of Celani et al. (Reference Celani, Cencini, Vergassola, Villermaux and Vincenzi2005). Second, inside the canopy layer (
$z< H$) the flow pattern can be very intricate with re-circulation regions behind the roughness elements and a channelling effect of strong flows along certain paths; this complex flow structure has a strong effect on scalar dispersion (Nosek et al. Reference Nosek, Kukacka, Kellnerova, Jurcakova and Janour2016; Castro et al. Reference Castro, Xie, Fuka, Robins, Carpentieri, Hayden, Hertwig and Coceal2017; Monnier et al. Reference Monnier, Goudarzi, Vinuesa and Wark2018). We expect that these flow patterns will become increasingly dominant in pair-dispersion statistics inside the canopy for separations of the order of the roughness obstacle widths 
$r \sim W$ (in our canopy 
$W=\frac {1}{2} H$). In our measurements, only the first approach towards this regime can be noticed in figure 5. Lastly, at scales much larger than the size of the canopy obstacles, 
$r\ggg H$, yet still inside the canopy, 
$z< H$, it could be expected that pair dispersion will be dominated by the random passage of the flow in between the roughness obstacles in an analogy to porous medium flow (Brenner Reference Brenner1980); this mechanical diffusion phenomenology (Nepf Reference Nepf1999) would lead to a dispersion regime with a pore scale equal to the obstacle spacing. Notably, for the time being these regimes remain a speculation intended to motivate future studies on this topic, as there are no previous pair-dispersion measurements in canopy flows at these scales.
Understanding pair dispersion is important for modelling dispersion in the environment, as it can be used to estimate the variance of passive scalar concentration fields from known sources. In particular, Cohen & Reynolds (Reference Cohen and Reynolds2000), proposed a Lagrangian stochastic modelling approach for pair dispersion in highly inhomogeneous turbulent flows such as canopy flows. Their approach is consistent with the necessary conditions introduced by Thomson (Reference Thomson1990), 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 (figure 8), we show explicitly at what time lag values this hypothesis is valid. In particular, 
$\rho _{v_{\| }}=0$ is reached for all tested initial separations at approximately 
$\tau \approx \tau _0$ (figure 8). This suggests that the concentration variance can be estimated from Cohen and Reynolds dispersion model for times 
$\tau > \tau _{v_{\| }} \approx \tau _b = ( {r_0^2}/{\epsilon } )^{1/3}$.
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 (Boffetta & Sokolov Reference Boffetta and Sokolov2002b; Bitane et al. Reference Bitane, Homann and Bec2012; Scatamacchia et al. Reference Scatamacchia, Biferale and Toschi2012; Shnapp & Liberzon Reference Shnapp and Liberzon2018; Tan & Ni Reference Tan and Ni2022), 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.
Acknowledgements
We are grateful to M. Konn, S. Kalenko, G. Gulitski, V. Babin, A. Shick and M. Hotovely, for their assistance in the wind-tunnel experiment.
Funding
This study is supported by the PAZY grants 2403170 and 1372020.
Declaration of interests
The authors report no conflict of interest.
Appendix. Subvolume flow parameters
Table 1. Turbulence parameters for each sub-volume for the 
$Re_\infty = 16\times 10^3$ case.
Table 2. Turbulence parameters for each sub-volume for the 
$Re_\infty = 26\times 10^3$ case.
Tables 1 and 2 present the flow parameters used in the calculations for each of the sub-volumes. Detailed calculations of these properties are given in Shnapp et al. (Reference Shnapp, Bohbot-Raviv, Liberzon and Fattal2020).
Open access
