Extending 3D-PTV for Lagrangian Measurements of Wind Tunnel Canopy Flows

We report on a unique extension to the three dimensional particle tracking velocimetry (3D-PTV) method. Its key feature is a simultaneous real-time particle segmentation image analysis during the high-speed digital video acquisition from multiple cameras. The real-time image analysis is managed on a dedicated hardware which outputs only the centroids, intensity and size of tracer particles. This results in a four to five orders of magnitude in data transfer rate reduction. As a result, the presented extension enables to work in large scale facilities and at high flow speeds, and to produce substantially longer experimental runs. We demonstrate the applicability of the proposed extension by presenting a unprecedented amount of Lagrangian trajectories recorded in an environmental wind tunnel experiment. Its importance is exemplified by revealing acceleration and dispersion statistics inside a heterogeneous canopy flow model through directly measured mean square displacements. These results are important for better understanding atmospheric turbulent flow dynamics and essential for constructing and validating Lagrangian dispersion models in the surface layer and especially in the roughness sublayer, where single point measurements are typically used along with the limitations of the Taylor “frozen turbulence” hypothesis.


Introduction
The Lagrangian representation is used to describe fluid properties along the positions of flow tracer particles, making it an important tool for analyzing turbulent transport and dispersion 1, 2 both in industrial and in environmental applications. Three dimensional particle tracking velocimetry (3D-PTV) is an experimental method which enables measuring and analyzing complex flows in the Lagrangian framework. 3D-PTV uses multi-view synchronized imaging to infer 3D positions of tracer particles carried by the flow 3 . These positions are linked in time and space (i.e., tracked) in order to construct trajectories -x(t) which represent the flow field in the Lagrangian framework. However, this well proven method can be hard to apply, partly due to vast amounts of data that need to be recorded, stored and handled.
High Reynolds number turbulent flows have an extremely broad dynamical range of length and time scales. Thus, a thorough sampling of a turbulent flow field has to be made at a high resolution to capture its finest details, and at the same time demands an extended sampling window to allow for a full description of all the scales of interest. This poses a practical limitation in the application of image based techniques, such as 3D-PTV, where high resolution images have to be acquired at sufficiently high rates and for prolonged durations of time, which result in a bottleneck for data transfer as well as data storage and analysis. This limitation is more severe in turbulent flows with strong mean components in which the acquisition rate threshold rises, and higher image resolution is required for tracking multiple tracers. Previous particle tracking methods were utilized to study homogeneous turbulence in a wind tunnel or water flume measurements [4][5][6] , turbulent pipe flow 7 , and turbulent boundary layers [8][9][10] . In some cases 4, 9 the camera(s) system can be installed on a traversing system moving at the mean flow speed, thus reducing to some extent the rate required during data acquisition. Nevertheless, traversing solutions are not applicable to most fluid flows that either exhibit strong shear or that are highly three-dimensional or inhomogeneous, such as in a mixing layer flow associated with two characteristic velocities. Real-time image compression techniques based on image binarization 11 , or edge detection 12 were suggested in the past to handle the bottleneck during data transfer. Unfortunately, these methods work only in ideal cases where tracers can be trivially segmented on a uniformly dark background. Such techniques do not work in experiments with realistic and more complex environments, where light reflections and unevenly illuminated background are inevitable. In this work we present an extension to 3D-PTV which is based on a new procedure that is different in hardware and software from the conventional 3D-PTV application. The essence of this extension is in transferring the image analysis step of the 3D-PTV algorithm to a real time hardware application, simultaneous with the data acquisition. As a result, a dramatic reduction is achieved in the required data transfer rate, as well as in storage volumes and in the computational time.
The new extension improves 3D-PTV for complex flows and enables converging Lagrangian statistics through significantly longer recording times.
Much of Earth's surface is accommodated with canopies of vegetation and buildings. These provide multi-scale sinks and sources of momentum and of scalars for the atmospheric boundary layer (in particular, water vapor, CO2 and heat). Understanding the role of forests in the global carbon cycle, monitoring pollution in high-raised mega-cities, as well as estimating momentum transfer within wind farms, are few examples associated with atmospheric turbulent canopy flows. Therefore, understanding the turbulent flow structure created by buildings and plants inside and above a canopy layer are important in many environmental aspects. They allow the development of closure schemes and parameterizations that account for surface layer contributions to the planetary boundary layer 13 and contribute in planning and optimizing field experiments and sampling.
The complexity of turbulent flows in canopies poses experimental challenges, both in the laboratory, and in the field. Laboratory experiments provide a well controlled environment for such studies with respect to environmental changes 14 (e.g. mean wind speed, direction and atmospheric stability). Such experiments are often conducted in wind tunnels or water flumes, in which obstacles of different shapes are used to model the canopy roughness. Measurements at fixed points are conducted using either hot wire, Laser Doppler Anemometry (LDA) [15][16][17][18][19][20][21][22][23][24][25][26][27] , or Particle image velocimetry (PIV) that provides details on the turbulent structure over two dimensional planes in the flow 8,[28][29][30][31][32][33][34][35] . Topological features can be studied with visualization techniques [36][37][38] , and with stereoscopic PIV 39 . All of the above measurement techniques are Eulerian based methods in which the flow is sampled at fixed points in space. Studies of canopy flows in the Lagrangian framework are rare. In a field experiment, DePaul & Shieh 40 have tracked neutrally buoyant balloons in an urban street canyon -highlighting a vortex cell parallel to the street canyon. More recently Di Bernardino et al. 8 have implemented particle tracking in a water flume above a two dimensional canopy model to obtain 2D Lagrangian temporal statistics. Inferring to canopy flow statistics from fixed point measurements requires a set of assumptions such as the Taylor "frozen turbulence" hypothesis, which are routinely violated in such non-homogeneous layers 41 . Furthermore, real canopies are three dimensional, and thus require measurements in 3D. 3D-PTV Lagrangian measurements can help to resolve these issues by providing multi-point measurements along three dimensional Lagrangian trajectories. Until now quantitative Lagrangian measurements within a 3D canopy layer have been unreachable, and have not been conducted in the past. Our first experimental attempt with limited success has been reported last year in a conference paper 42 .
The extension we propose here to the 3D-PTV method is applied to investigate a canopy flow model within a full-scale environmental wind tunnel. This implementation required laser illumination in the forward scattering mode, systematic calibration using a remotely operated robotic arm, and an elaborated post-processing using the open source openPTV 43 software and the Flowtracks 44 package. This combined setup enabled to successfully reconstruct millions of trajectories within and above the canopy flow model. Using this dataset we explore previously unattainable Lagrangian statistics -in particular three dimensional Lagrangian accelerations and dispersion within the canopy. Our initial analysis reveals that the mean shear above the canopy 45,46,48 affects the standard deviation of Lagrangian acceleration. Furthermore, we find an analogy of the variance of displacements with the dispersion theory of Taylor 47 . The diffusivity that we find is governed by the standard deviation of velocity and a length scale that is independent of space and seemingly related to the turbulence in the wake of roughness elements.
inside the canopy the trajectories are highly curved, and are moving in all directions. This visualization highlights the strong inhomogeneity of the flow in the canopy roughness layer.

Velocity Distributions
Probability distribution functions (PDFs) of the instantaneous streamwise velocity component for the Re ∞ = 2.6 × 10 4 case, normalized with U ∞ , are presented in Fig. 2. Values for the PDFs are sampled from 5 different heights over an entire canopy roughness cell, thus include phenomena related to spatial variation of the velocity field. Therefore, the PDFs are not directly comparable with the more common Eulerian single fixed point PDFs of the turbulent velocity u . The modes (most probable values) and the widths of the PDFs increase while rising from within the canopy layer to the roughness sublayer right above the canopy z ≥ H, covering velocities higher than the free stream wind velocity, U ∞ , and demonstrating positive skewness. Positive skewness of the streamwise velocity is characteristic of canopy flows 48 , typically associated with the asymmetry of high velocity sweeps and low velocity ejections 49,50 . However, in Fig. 2 the positive skewness could arise for different reasonspartially due to the spatial averaging, and partly to a known bias in imaging based techniques that over-sample low velocity tracers. This over-sampling is related to slow particles that remain longer in the field of view, providing longer trajectories, and adding weight to samples of lower velocities 3 .In addition, the seeding nozzles were not equally distributed in the wind tunnel cross-section, being closer to the bottom floor; this might have resulted in the flow tracers preferentially sampling low momentum flow from the bottom wall as compared to the high velocity of the free stream. Our future study will be devoted to a careful comparison of these effects using several measurement methods such as particle image velocimetry and Laser Doppler velocimetry, in addition to 3D-PTV, in the same canopy setup.

Lagrangian Accelerations
The Lagrangian acceleration, a = D v/Dt, is a small scale turbulent property 51 , characterized by a short correlation time, on the order of the Kolmogorov time scale 10,51,52 . Due to the separation of scales in high Reynolds number turbulent flows, it was assumed that acceleration statistics are independent of the large scale structures of the flow. We address this hypothesis in the case of the canopy flow model.
We estimate statistics separately for the three components of Lagrangian acceleration, sampled at 450 locations. At each location we use trajectories that cross a sphere of radius 0.05H (the size of the sample volumes and number of sample locations have do not affect the results, these validation steps are not shown for the sake of brevity, each sphere contains at least 1000 samples). The PDFs of acceleration components are constructed from the results of all the sampling locations, utilizing flow stationarity and ergodicity 2 . In Fig. 3 we present PDFs of the three components of Lagrangian accelerations, normalized by the respective standard deviations, for the the two Re ∞ cases. The scattered data from random sample locations is added to the plot to represent the width of the distribution. The acceleration components PDFs exhibit long tails characteristic of the intermittent nature of the Lagrangian accelerations [52][53][54] . The data for all the studied cases collapse on the same curve and statistics are remarkably similar regardless the inhomogeneity of the flow in the canopy. Furthermore, the stretched exponential function The solid curve corresponds to Eq. (1) with parameters β = 0.45, σ = 0.85, γ = 1.6 and C = 0.525. The obtained exponent, γ, agrees with the zero mean shear flow cases 54,55 . This fact suggests that the normalized Lagrangian acceleration PDFs are insensitive to the inhomogeneity and strong mean shear at the canopy roughness sublayer. Following the standard approach to canopy flows analysis 48 , we obtain horizontally averaged statistics as a function of height. In Fig. 4 we present the horizontally averaged acceleration standard deviationσ a i -as a function of height for the three acceleration components and the two Re ∞ cases. The range of uncertainties, presented as error bars, is estimated using the bootstrapping method where the data is divided into subsets and error propagation is integrated in the horizontal averaging. For all the profiles, σ a i decreases with the decrease in height, from 1.5H toward the canopy height down to about 1.25H, where there is a local increase with a notable maximum at z ≈ H, and further decreases inside the canopy z < H. Notably, the local maximum of Lagrangian acceleration variance profile coincides with previously reported peak in turbulent kinetic energy dissipation profiles 56 at the height typically associated with the canopy roughness sublayer. The shape of the profiles is similar for the two Re ∞ cases, yet does not collapse on a single curve with a straightforward normalization.  In homogeneous isotropic turbulence Kolmogorov similarity theory leads to a scaling of the variance of Lagrangian accelerations 2σ 2 a i = a 0 ε 3/2 ν −1/2 , where ε is the mean turbulent kinetic energy dissipation rate, ν is the kinematic viscosity, and a 0 is a Re dependent coefficient 51,54,57,58 . Using ε = σ 3 u /L with σ 2 u = (σ 2 u x + σ 2 u y + σ 2 u z )/3, and L the integral length scale 59 , one arrives at: Previous works in zero mean shear turbulent flows presented σ a i versus σ u 54 , or conditional acceleration variance a 2 |σ u 60 and found good agreement with the 9/2 power law in Eq. (2). In Fig. 5 we plot the streamwise acceleration variance σ 2 a x against σ 9/2 u for the two Re ∞ cases. Similarly to the zero mean shear case 54, 60 , we observe that σ a x scales with σ where σ 2 a x decreases with σ u locally. The same behavior is seen in both Re ∞ cases, and were observed for all acceleration components (not shown for brevity). The observation of two different regions where σ 2 a x ∝ σ 9/2 u are in line with the proposed "mixing layer" structure of the canopy flows 20,45,48 . The two trends in Fig. 5 suggest that the two regions of two different turbulent flows interact through a a shear layer interface at the top of the canopy, H ≤ z ≤ 1.25H. We suggest that according to Eq. (2) the two slopes provide evidence that the small scale structure of the turbulence is different inside and outside of the canopy layer, due to different magnitudes of the mean shear, and the different degrees of inhomogeneity in these regions.

Dispersion
In an inhomogeneous canopy flow field, the velocity variance and the auto-covariance required for dispersion estimate according to Taylor theory 59 vary in space. However, the unique dataset enables us to estimate turbulent dispersion in a canopy flow directly using the Lagrangian trajectories. We examine Lagrangian tracers in respect to their first observed positions, x 0 (t 0 ), and note that for stationary flows the statistics depend on time difference τ = t − t 0 2 . Furthermore, we define the ensemble average · over all trajectories passing through x 0 , with respect to the same time difference τ. Thus we define the mean particle displacement as: The transverse dispersion is obtained as the variance in respect to the mean displacement ∆x i ( x 0 , τ): Finally, we define the diffusivity K as one half the time derivative of the dispersion 59 : We present results for three positions x 0 , located on a single vertical line at three different heights for the case of Re ∞ = 1.6 × 10 4 . In Fig. 6(a) we present the displacements, ∆x(τ) = x(τ) − x 0 , of 1000 randomly chosen trajectories (think black lines) and also the ensemble average displacement in the streamwise direction. The dispersion is reflected in the width covered by trajectories relative to the ensemble average displacement curve (thick dashed line); while the streamwise velocity of trajectories is proportional to the slope of each line. The finite size of the measurement volume limits the range that we could estimate for the dispersion. High free stream flow velocity (for example the highest point seen in Fig. 6), leads to faster decreasing variance as compared to the results obtained at lower positions and lower velocities.   4), at the three locations as a function of τ. The dispersion is largest at the z = 1.5H and decreases with the height. Interestingly, for all the three heights we observe a similar behavior of an initial quadratic growth in time (reflecting a ballistic regime) followed by a linear growth, suggesting a constant diffusivity, Eq. (5). The slopes decrease slightly at the large τ, due to particles leaving the measurement volume.
In the inset of Fig. 6(b) we plot the dispersion from three different heights that collapse on a single curve if we scale the time differences with the standard deviation of streamwise velocity at the appropriate heightσ u x ( x 0 ). This result is robust for measurements at other points and other Re ∞ case (not shown for the sake of brevity). Using Eq. (5) one can express the diffusivity in terms of a length scale L , K T,x = 1 2 L σ u x , noting that the slope in the inset of Fig. 6(b) is equal to 1 2 L . The best fit shown in the inset gives L = 4.1 ± 0.1 mm. This observation can be seen as an analogy to previously reported relation between turbulent kinetic energy and diffusivity in the two dimensional turbulence 61 , and with the constant "mixing length model" used for momentum transfer in boundary layers 46 . This observation can also be interpreted as a constant drag length scale 62 , possibly related to turbulent fluctuations in the wake of the canopy obstacles which thickness is close to L . It is also remarkable that the dissipation, shown in Fig. 5 changes with the height, unlike the length scale that is apparently independent of height, as seen in Fig. 6(b).

Discussion
We report on an extension to the 3D-PTV method based on an extensive real time image analysis on a dedicated hardware and software. This extended 3D-PTV scheme has several advantages: (i) It makes it possible to conduct very long experimental runs under difficult imaging conditions, which are demonstrated here on a canopy flow model in a full scale environmental wind tunnel. (ii) It reduces the time required for image analysis, thus increasing by several orders of magnitude the amount of recorded data that is valuable for the post-analysis. (iii) It allows to preform three-dimensional measurements in turbulent flows where high spatial and temporal resolutions are essential, without compromising on the duration of experimental runs.
The proposed method has been successfully applied to measure for the first time Lagrangian trajectories within a canopy layer modeled in a full scale environmental wind tunnel. This opens up an extraordinary opportunity to explore the inhomogeneous roughness sub-layer in the Lagrangian framework. The intermittent structure of turbulence in the canopy layer requires large 7/15 data set to obtain converged statistics which is impossible to obtain from short recording sessions.
We obtained unique Lagrangian turbulent flow statistics in the canopy flow model. We found that the standardized Lagrangian acceleration PDFs shape is independent of position, direction, and the free stream velocity in the range tested, and have a shape similar to the previously measured homogeneous isotropic turbulent flow cases. The amplitude of the acceleration, described here with the standard deviation of Lagrangian accelerations, increases with height above the bottom wall, and has a local maximum at z ≈ H. More interestingly, the relation between Lagrangian acceleration standard deviation to the root mean square of turbulent velocity seems to agree with Kolmogorov scaling but in two distinct regions, inside and above the canopy. Between the two different scalings we observe a transition region in which turbulent kinetic energy is not related to the local (in space) dissipation and small scales, in agreement with the mixing layer analogy suggested for the canopy flows 20,45 . Furthermore, the Lagrangian data also enables to examine dispersion through the direct measurements of particle displacements passing through spatial locations at different heights. Our analysis reveals that there is a ballistic type of dispersion followed by a constant diffusivity range. The diffusivity, defined in a form analogous to Taylor's 47 dispersion theory, scales with the turbulent velocity standard deviation. Furthermore, it can be presented in a form of a length scale that seems to be predefined by the canopy obstacles and not by the mean shear.
The measurements presented here hold a lot of additional information that need to be processed and analyzed. The present results are only the first attempt to highlight the possibilities opened by the developed 3D-PTV extension, and the importance of three-dimensional Lagrangian studies in a turbulent canopy flow.

3D-PTV extension
The infographic in Fig. 7 presents the main concepts of the 3D-PTV extension. It draws the line through the history of 3D-PTV, highlighting the past, present and future development milestones. In traditional 3D-PTV analog or digital video is streamed from an imaging device to a storage media at a maximum possible transfer rate, which is typically limited by the recording speed of a storage media. The first generation of 3D-PTV systems were based on 640 × 480 pixels videos (VGA format) 3 , whereas the present high speed digital video reach 4 million pixel frames at a recording rate of 500 frames per second. The present maximal data transfer rates are about 2.4 Gb/s if video is streamed to a digital video recording system, writing the stream in parallel into multiple hard drives using high throughput cables and dedicated electronics. As shown in red boxes in Fig.fig:infographic, if a recording experiment lasts 1 hour using a 4-camera 3D-PTV system, a digital recording system will transfer and store about 7 × 10 6 images equivalent to about 200 Terabytes. Post processing of these images into Lagrangian trajectories can then take an estimated period of 5 months, including the image reconstruction from a fast binary format, image segmentation analysis (mainly detection of objects, assumed O(1) sec per frame), stereo-matching and Lagrangian tracking. In 2007 and 2011, simple image processing algorithms pioneered the compression of transferred images 11,12 (blue boxes in Fig.fig:infographic). With image compression, the data from the 1 hour experiment can be compressed to a total of ≈ 200 Gb. Although it significantly reduces the amount of recorded data, it does not resolve the time-consuming post-processing steps that would take about a 1 month of work. Furthermore, a main limitation of these simple image compression algorithms is that they cannot perform well in complex image environments with light reflections and uneven illumination in the background. Unlike image compression, in the proposed extension, the system performs complex image analysis (presented in details below and in the supplementary materials) and as a result streams only 2D blob-coordinates. This enables direct writing onto an off-the-shelf hard drive that can store the data at such a rate for many days. For instance, the 1 hour recording (equivalent to 200 Tb of video data) is processed and transformed into a ∼ 20Gb data set, like in the presented 3D-PTV canopy flow experiment. The real time image analysis extension yields a 1 : 10000 compression ratio of the data which turns out to be essential for studying canopy flow models and other wind tunnel experiments which up to date could not have been reached by 3D-PTV. The unique 3D-PTV system lends access to Lagrangian information of the flow in the inter-obstacle regions through flexible positioning of the measurement volumes. It can also record data for hours, quantifying the changes in the flow both on very short time scales and for many turnover time scales. Independently of the length of recording and number of runs, the system presents an important advantage of the ability to post-process the data and obtain Lagrangian trajectories during the same day. It allows experimentalists to improve and optimize the settings for the next day, increasing the repeatability of the experiments. Fig. 7 also presents the future 3D-PTV generation that allows complete 4D, real-time measurement technique. In principle, such real-time 4D-PTV could be achieved by formulating/operating real-time stereo-matching and Lagrangian tracking algorithms, presumably on the same hardware.

Real-time image analysis on hardware
The proposed extension is designed to resolve the data transfer bottleneck using the real-time image analysis system (called internally "blob recorded") based on high-end framegrabbers (microEnable v5) with FPGA hardware and software (Silicon Software GmbH, Germany). The system was designed and integrated in collaboration with one of the authors E.S., a specialists 8/15 Figure 7. Infographic diagram presenting the development of 3D-PTV imaging and recording systems along three generations. During this period, the advance in imaging technology enabled to grow from VGA format images at 30 frames per second, to four million pixel of recorded images at 500 frames per second, limited by the 2.4 Gb/s data transfer rate. In the straightforward approach (red boxes), images are directly stored on a hard disk and analyzed in the next step resulting in vast amounts of data that have to be stored and managed. Real-time image compression 11, 12 (blue boxes) enabled to reduce this large amount of stored data, as well as the computational time under restricted idealized imaging conditions. The elaborated real-time image analysis, introduced in this work (green boxes), enables a substantial reduction in computing times and storage volumes, and has been developed to perform equally well under non-ideal imaging conditions with strong reflections and uneven illumination, such as in wind tunnels. The numbers indicate estimated storage volume and computational time for each PTV generation, starting from a 1 hour recording experiment. The suggested extension opens the gate for full real time tracking that might be formulated on FPGA cards in the near future (orange box). from 1Vision LTD (Netanya, Israel). The system runs complex image analysis on acquired images in real time and at a high frame rate. The output is a data stream of the position from each detected object in pixel coordinates, along with its size and a bounding box. Real-time image analysis provides the 3D-PTV method with a real-time input simultaneously with the image acquisition, a step that was traditionally performed off-line on stored images. Subsequently, the tracer's positions in pixel units are incorporated with the 3D camera calibration, stereo-matching and tracking in order to facilitate 3D Lagrangian trajectories 3 . A block diagram of the algorithm used for real-time image analysis, is shown in Fig. 8. It is basically a customized object identification and analyzing method, commonly implemented in image processing softwares, that could only recently be performed in real-time on a FPGA with a sufficient memory and computational capabilities. A detailed description of the blob-recorder can be found in the Supplementary Materials.

The wind tunnel canopy model
We preformed the experiment in the Environmental Wind Tunnel Laboratory at the Israel Institute for Biological Research (IIBR). The facility is an open circuit wind tunnel featuring a 14m × 2 × 2 m 2 test section 42 . A surface layer typical of a dense canopy flow was created by placing "L" shaped rectangular roughness elements along the tunnel's floor, and four spires 15 upstream to the measuring section, as shown schematically in Fig. 9. As usual, we set x to coincide with the streamwise direction, y as horizontal cross-stream and z as the the wall normal direction pointing against gravity. We introduced heterogeneity in the canopy by using roughness elements of two heights -0.5H and 1H, where H = 100 mm. The element's width was 0.5H and the thickness was 4 mm. We arranged the elements in a staggered configuration. Each lateral-row of elements was kept uniform in height, whereas alternating the height between consecutive rows. In the first few meters of the test section elements were fitted at a low density with spacing of 2H. Starting 50H upstream from the measurements location, a dens, staggered, double-height canopy layer was placed, with lateral spacing of 0.5H and a streamwise distance between consecutive rows of 0.75H. The canopy frontal area density, defined as λ f = A f /A T , (where A f is the frontal area of the elements and A T is the lot area of the canopy layer), was 9 16 . The experiments were conducted at nominal free wind speeds of U ∞ = 2.5 and 4 m s −1 ,

Seeding
Hollow glass micro-spheres (Potters Industries, Sphericell) were used as tracers. A scanning electron microscope inspection confirmed that the tracers are spherical with diameters in the range [2 − 25] µm with a mean diameter of d p ≈ 10 µm. The mean particle response time in the Stokes regime is estimated as τ p = ρ p d 2 p /18µ ≈ 0.4 ms, based on d p and ρ p = 1000 kg m −3 . Following previous studies 63 , the highest flow frequency resolved by the particles is 500 Hz (to an error of 5% in velocity RMS). Based on an assessment of the Kolmogorov time scale ∼ 7 ms, this value should be sufficient for the particles to resolve the turbulence dynamics, yet a minor filtering effect of the highest frequency content is possible. The particles were seeded in the flow from four point sources located at two different heights, 7H and 3H and placed 60H upstream from the measurement position (see Fig 9). To reduce agglomerations due to humidity, the tracers were dried and stored in an oven for 24 hours prior to seeding. A great deal of effort was invested in constructing the pneumatic seeding devices and on their positioning and air-intensity, to ensure a well mixed and uniform seeding at the measurement location.

Experimental setup
Intense illumination was required to identify the fast moving small tracer particles in the images. For this reason we used a 10 Watt 532 nm continuous wave laser (CNI lasers, MGL-V-532). The laser beam was expanded to an ellipsoidal shape with radii 80 × 40 mm through a pair of cylindrical lenses, rotated 90 degrees relative to each other around the beam axis. The laser was positioned outside the test section on the one side of the wind tunnel. Four high-speed CMOS cameras (Optronis CP80-4-M/C-500, 100 mm lenses) were positioned on the opposite side of the test section, pointing towards the measurement location (see Fig. 9). Thus the tracers were illuminated in the forward scattering mode that is stronger than the more commonly used side scatter mode. The strong illumination enabled to reduce the shutter time and improve the focal depth of the cameras using larger F-number lenses. The increased number of out-of-focus tracers were filtered out by the real-time image processing through an adaptive threshold algorithm.
The volume observed by the cameras at the measurement location and the actual spatial resolution were approximately 80 × 40 × 40 mm, (in streamwise, wall normal and span-wise directions) and 40µm per pixel, respectively. The measurements were conducted at 20 sub-volumes distributed over 5 heights in order to sample the complete single canopy unit-cell within the staggered, double-height canopy layer.
Calibration is an essential part of the 3D-PTV method that dictates, to a large extent, its uncertainty in the determination of the tracer positions 3 . In the wind tunnel we used a three-dimensional calibration target composed of 80 white dots, each of 0.5 mm in diameter, at known positions over 3 vertical planes. The target was mounted on a remotely controlled three-axis traverse system with a positioning accuracy less than 40 µm. The dense canopy layer prevented us from accessing the measurement location during the period of measurements. Thus, experimental repeatability was achieved through the use of the high precision traverse system allowing multiple calibrations prior to and after each experimental run. Typical calibration errors in this experiment are estimated based on repeated calibration and detection of the calibration targets as 15 µm in the x, z (streamwisefloor-normal, see Fig. 9) plane, and 30 µm in the y, along the imaging axis (depth). At the highest frame rate of 1000 Hz, this error could propagate into a single point velocity uncertainty of approximately 15 mm s −1 , and 30 mm s −1 , respectively. Note that generally in a 3D-PTV experiment the velocity uncertainty is estimated using additional available information in space and time, along the Lagrangian trajectories 4 .
The recording rate was 500 Hz at the full resolution of 2304 × 1720 pixels inside the canopy surface layer, and 1000 Hz at a reduced resolution 2304 × 860 pixels above the canopy. The settings are a trade-off between the requirements for illumination, higher recording rate due to large velocities and shorter exposure times (to avoid smearing of tracer images). The cameras were synchronized with an external timing controller (CC320, Gardasoft). The maximal error in timing is significantly shorter than a 1/frame rate (∼ 10 µs). The recording system is responsible for counting the recorded frames as a sanity check, and each frame is recorded with a unique time stamp.
The data was gathered in 40 individual experimental runs, each lasting between 10 to 15 minutes. These time intervals, being roughly in the order of 2 × 10 3 turnover time, ensure sufficient sampling of statistically independent trajectories. The turnover time is estimated as H/ u H , where for the high flow rate, u H = 0.28 m s −1 is the mean (spatial and temporal) streamwise velocity at z = H.

Post Processing
We used the predictor-corrector 4 frame steps tracking algorithm of OpenPTV 64 , that was successfully used in various turbulent flow applications 4, 65-68 .Due to the flow inhomogeneity we followed a systematic routine for choosing the different tracking parameters at the different measurement sub-volumes, similar to the previous work 67 . In each region of the flow, we chose a sub-sample of the data and applied tracking with different tracking parameters. The parameters are tuned until the estimated root mean square of velocity, the number of trajectories, and the length of trajectories reached a plateau. The software modifications are available from the open source repository, OpenPTV 43 .
The use of real time imaging prevents the usage of the space-image tracking algorithm 69 , that requires access to the original images. Instead, we post-process the trajectories using the position-velocity space linking algorithm 70 . The trajectories data was filtered and differentiated using the standard method of OpenPTV 44,66 . This consists of fitting a second order polynomial to a window of N frames for each three dimensional trajectory. The polynomial is used also to estimate first and second order time derivatives of the positions, i.e velocity and acceleration. The data at the edges of the trajectories, where the N sized window cannot be fitted symmetrically, were discarded. The cutoff frequency of this low pass filter is f c = 150 Hz. Post processing was conducted using the OpenPTV post-processing package, called Flowtracks 44 . It includes operations such as storing, handling and manipulating capabilities of the vast amount of data through a sophisticated application of the HDF5 format.

Data availability
The processed datasets from the current study could be obtained from the corresponding author.