Turbulence -- Obstacle Interactions in the Lagrangian Framework: Applications for Stochastic Modeling in Canopy Flows
Ron Shnapp, Yardena Bohbot-Raviv, Alex Liberzon, Eyal Fattal
TTurbulence – Obstacle Interactions in the Lagrangian Framework: Applicationsfor Stochastic Modeling in Canopy Flows ∗ Ron Shnapp, † Yardena Bohbot-Raviv, ‡ Alex Liberzon, † and Eyal Fattal ‡ Lagrangian stochastic models are widely used to predict and analyze turbulent dispersion in com-plex environments, such as in various terrestrial and marine canopy flows. However, due to a lackof empirical data, it is still not understood how particular features of highly inhomogeneous canopyflows affect the Lagrangian statistics. In this work, we study Lagrangian short time statistics byanalyzing empirical Lagrangian trajectories in sub-volumes of space that are small in comparisonwith the canopy height. For the analysis we used 3D Lagrangian trajectories measured in a densecanopy flow model in a wind-tunnel, using an extended version of real-time 3D particle trackingvelocimetry (3D-PTV). One of our key results is that the random turbulent fluctuations due tothe intense dissipation were more dominant than the flow’s inhomogeneity in affecting the short-time Lagrangian statistics. This amounts to a so-called quasi-homogeneous regime of Lagrangianstatistics at small scales. Using the Lagrangian dataset we calculate the Lagrangian autocorrelationfunction and the second-order Lagrangian structure-function, and extract associated parameters,namely a Lagrangian velocity decorrelation timescale, T i , and the Kolmogorov constant, C . Wedemonstrate that in the quasi-homogeneous regime, both these functions are well represented usinga second-order Lagrangian stochastic model that was designed for homogeneous flows. Further-more, we show that the spatial variations of the Lagrangian separation of scales, T i /τ η , and theKolmogorov constant, C , cannot be explained by the variation of the Reynolds number, Re λ , inspace, and that T i /τ η was small as compared with homogeneous turbulence predictions at similar Re λ . We thus hypothesize that these characteristics occurred due to the injection of kinetic energyat small scales due to the so-called “wake production” process, and show empirical results support-ing our hypothesis. These findings shed light on key features of Lagrangian statistics in flows withintense dissipation, and have direct implications for modeling short term dispersion in such complexenvironments. I. INTRODUCTION
Scalar dispersion in the atmospheric surface layer is strongly influenced by the turbulent canopy flows.These flows affect the dispersal of pathogens and the ventilation of urban areas [1], as well as the dispersal ofspores, bacteria, and seeds in forests and fields [2]. After years of research, it had become commonly acceptedthat there is a particular difficulty to model dispersion in canopy flows through Eulerian models, resultingfrom the failure of gradient-diffusion theory in these flows [3–5]. Consequently, Lagrangian Stochastic Models(LSM) gained popularity as a state-of-the-art modeling approach within the community (e.g. Refs. [2, 6–16]among others). LSMs can efficiently predict transport and dispersion in turbulent flows through Monte-Carlosimulations [8], and are specifically useful in applications to inhomogeneous turbulence (i.e. where the flowstatistics vary in space) and to cases with complex distributions of sources.In the LSM framework, a stochastic Markov random walk provides increments of the states of “marked”particles, defined by the particles’ position and its time derivatives up to n order d n x /dt n , (usually n = 1 or n = 2 is the order of the model, and boldfaced symbols denote vectors). There are two fundamental steps inconstructing an LSM: a) determine an appropriate stochastic process defining the random walk, and b) toobtain functional expressions relating the Lagrangian operators of the random walk to statistics of the flowfield. In other words, one must find Lagrangian equations of motion in terms of a priory known Eulerianvelocity statistics.There are significant challenges associated with the two steps above due to the complex nature of thecanopy flows, in which various processes occur simultaneously at different scales. In particular, there is anumber of topics that demand special attention when constructing LSMs for canopy flows:1. turbulence generated simultaneously by both the large scale shear and due to flow–obstacle interac-tions (wakes) by canopy elements, leading to the so-called spectral short-circuiting through wake– ∗ [email protected] † School of Mechanical Engineering, Tel Aviv University, Tel Aviv, Israel ‡ Israel Institute for Biological Research, Ness Ziona, Israel a r X i v : . [ phy s i c s . f l u - dyn ] J u l production [17, 18];2. non uniqueness in the formulations of LSMs for inhomogeneous flows, namely in the solution of Thom-son’s well-mixed principle [19, 20];3. inhomogeneity due to the significant difference between the flow above and inside the canopy layer [15],leading to a strong shear layer near the canopy top;4. Uncertainty in the parameterizations of Lagrangian statistics that are needed for LSM construction,for example the Kolmogorof constant or the Lagrangian integral timescale [21–23];5. the prevalence of large coherent structures due to the phenomenon so-called the mixing layer anal-ogy [24], leading to non-Gaussian distributions of the Eulerian velocity [25];6. the “mechanical diffusion” as a result of the fluid having to bypass the canopy obstacles [26] (similarto dispersion in porous media flows);Notably, the references above could not have been a conclusive list due to the vast body of literature. Inthis work, we focus on small scale motion in wakes of canopy obstacles associated with item 1. Specifically,turbulent kinetic energy is produced in canopies in two ways: production due to mean shear at scalescomparable with obstacle height ( H ), and production due to obstacle drag at smaller scales, so-called wake-production. Wake-production is said to short-circuit the turbulent cascade since it extracts energy from largescales and injects it directly at smaller scales [17]. However, except for the theory by Poggi et al. [18] showingthat wake production can affect the Kolmogorov constant C (introduced below), how wake productionaffects Lagrangian statistics in canopy flows is not known. Importantly, direct estimations of Lagrangianflow statistics in canopy flows that could provide crucial data to resolve the items above are lacking.Technological and scientific advances of the last two decades enabled the gathering of invaluable empiricaldata in the Lagrangian framework, both through experiments and direct numerical simulations (e.g. Refs. [27–46]). Such studies have analyzed Lagrangian dynamics in details, bringing attention in the community todelicate mechanisms underlying the motion of Lagrangian particles in turbulent flows [40]. The empiricaldata obtained in such studies can be used to bridge the gap between the fluid mechanics and the stochasticmodels through more accurate parameterizations and validation of theories. In most cases, due to practicallimitations, the focus was put on homogeneous isotropic turbulent flows (HIT); however, in the recent years,studies with embedded inhomogeneity and anisotropy are gaining more attention (e.g. Refs. [43–45, 47–49]).In this study, we use measurements of flow tracers’ trajectories in a heterogeneous canopy modeled in awind tunnel, [49, 50], to better understand the effect of flow–obstacle interactions on short time Lagrangianstatistics (item 1). We focus our analysis on small scales by exploring Lagrangian statistics inside small sub-volumes of space, namely smaller than (cid:104) u (cid:105) ∂ (cid:104) u (cid:105) /∂z where (cid:104) u (cid:105) is the streamwise mean velocity. Our analysis revealsthat the short time Lagrangian statistics were affected more strongly by random turbulent fluctuations than bythe flow inhomogeneity due to intense dissipation; this amounts to quasi-homogeneous regime of Lagrangiandynamics at short-time and small-scale. Accordingly, the Lagrangian autocorrelation function and the second-order structure-function could be well represented by a second order LSM designed for homogeneous flows.Furthermore, our empirical study in the quasi-homogeneous regime suggests that Lagrangian separationof scales and the Kolmogorov constant (definitions below) are affected by wake-production, in agreementwith Poggi et al. [18]. The results of our study are relevant for flows with intense dissipation and spectralshort-circuiting effects. Importantly, we demonstrate how measurements performed directly in the Lagrangianframework can provide crucial insight into Lagrangian dynamics in canopies, and thus that they have potentialfor improving dispersion models’ accuracy in canopy flows.The rest of this paper is organized as follows. In Section II we present three LSMs that will be usedin our analysis, and the details of our experiment and our analysis. In Section III, we reveal a quasi-homogeneous regime of Lagrangian statistics in the small scales. In Section IV we present direct estimationsof the Lagrangian velocity autocorrelation, the structure-functions and the associated parameters ( T i and C ),demonstrating their spatial distributions. In Section V, we compare our empirical results with predictionsfor homogeneous flows and hypothesize on the role of wake production. Lastly, we summarize and presentour conclusions in Section VI. II. METHODSA. Definitions
Let us first define three different LSMs that will be used in the consequent analysis. The simplest first orderLSM ( n = 1) assumes that the turbulence has spatially homogeneous statistics and that it is characterizedby Gaussian velocity PDFs at each point (hereafter called Gaussian turbulence). This model is essentiallythe Ornstein-Uhlenbeck (OU) process [8, 51]: dv (cid:48) i = α i ( x , v , t ) dt + β ij ( x , v , t ) dξ j ; dx i = [ v (cid:48) i + (cid:104) u i (cid:105) ] dτ (1) α i = − v (cid:48) i T L,i ; β ij = δ ij (cid:18) σ i T L,i (cid:19) / . (2)Here v i = (cid:104) u i (cid:105) + v (cid:48) i is the i th component of the velocity of a Lagrangian particle where (cid:104) u i (cid:105) is the meanvelocity field and v (cid:48) i a velocity fluctuation, dξ i are increment components of a Wiener process, δ ij is theKronecker delta function, and σ i is the standard deviation of v (cid:48) i . The Lagrangian velocity autocorrelationfunction, ρ ij , and the Lagrangian integral timescale, T L,i , are defined as: ρ ij ( τ ) = (cid:104) v (cid:48) i (0) v (cid:48) j ( τ ) (cid:105)(cid:104) v (cid:48) i (0) (cid:105) / (cid:104) v (cid:48) j ( τ ) (cid:105) / ; T L,i = (cid:90) ∞ ρ ii ( τ ) dτ (3)where angular brackets (cid:104)·(cid:105) denote an average, the velocities v i (0) v j ( τ ) are taken along an individual La-grangian trajectory, and since we are considering a homogeneous and stationary flow, ρ ij ( τ ) is a functionof τ only. Doob [52] proved that the OU process, Eq. (1), is essentially the only process defined with theproperties: stationary, Gaussian and Markovian, with an exponential autocorrelation function: ρ ij ( τ ) = δ ij exp (cid:18) − τT L,i (cid:19) . (4)Furthermore, based on the Obukhov conjecture [53], consistency with the Kolmogorov inertial range scalingof Lagrangian velocity increments [54] requires that [8] –2 σ i T L,i = D ii ( τ ) τ = C (cid:15) ; D ii ( τ ) = (cid:104) [ v (cid:48) i ( τ ) − v (cid:48) i (0)] (cid:105) (5)where D ii ( τ ) is termed the Lagrangian second order structure function, C is the so-called Kolmogorovconstant, and (cid:15) is the mean rate of turbulent kinetic energy dissipation. Borgas and Sawford [55] showedthat for HIT these choices of α i , β i are unique solutions of Thomson’s well-mixed condition [19] and thus itis an exact solution of the Fokker-Plank equation.In realistic situations the flows’ statistics vary in space, namely flows are inhomogeneous, and this fact isnot treated by the simplest model above, eq. (2). To overcome this issue, Thomson [19] solved the well-mixedcondition and provided a LSM for such flows assuming that the turbulence is Gaussian and consistent withKolmogorov similarity for D ij . This Markovian LSM for inhomogeneous flows still uses eq. (1), but thecoefficients change according to: α i = − C (cid:15) R − ij v (cid:48) j + φ i g ; b ij = (cid:112) C (cid:15)δ ij φ i g = 12 ∂R il ∂x l + 12 R − lj (cid:104) u k (cid:105) ∂R il ∂x k v (cid:48) j + 12 R − lj ∂R il ∂x k v (cid:48) j v (cid:48) k . (6)Following Wilson et al. [56], this equation can be written more compactly as: φ i g = T i + T ij v (cid:48) j + T ijk v (cid:48) j v (cid:48) k (7)where the T m are coefficients multiplying the fluctuating velocity to power m . As written in the introduction(item 2), this model is not a unique solution of the well-mixed condition under the above assumptions.Nevertheless, it is commonly used for modeling dispersion in canopy flows (e.g., Refs [10, 23, 56, 57]) and weshall utilize it here as well.So far we have only considered first order LSMs ( n = 1) that are regarded accurate in flows with veryhigh Reynolds numbers. In such cases the particle accelerations have very short correlation times relativeto the integral timescales [27, 28, 54]. However, if this separation of timescales reduces, for example due tofinite Reynolds number effects, the first order models become less accurate in modeling the dispersion. Theinfluence of finite separation of scales was addressed by Sawford [58] in the introduction of a second ordermodel ( n = 2) for a homogeneous turbulence case: T L da i + (1 + R / ) a i ( t ) dt + R / T L (cid:90) t a i ( τ ) dτ dt = (cid:114) σT L R / (cid:0) R − / (cid:1) dξ i dv i = a i dτ ; dx i = ( v i + (cid:104) u i (cid:105) ) dτ ; R = 16 a C (cid:18) τ e τ η (cid:19) (8)where a is a Kolmogorov constant for the variance of acceleration, a i is a component of the particle’sacceleration vector, τ e is an Eulerian integral timescale, and τ η = ( ν/(cid:15) ) / is the Kolmogorov timescale. Themodel Eq. (8) has the following velocity autocorrelation function [58]: ρ ij ( τ ) = δ ij T ,i exp( − τ /T ,i ) − T ,i exp( − τ /T ,i ) T ,i − T ,i ; T L,i = T ,i + T ,i . (9)This second order model was extended by Du et al. [59] to the decaying turbulence case and by Reynolds [20]to include vertical inhomogeneity of the turbulent flow statistics.Above we have considered three models: Eq. (1) with Eq. (2) which is a first order model for infiniteReynolds number homogeneous turbulence, Eq. (1) with Eq. (7) which is also a first order model but forinhomogeneous flows, and lastly, Eq. (8), which is a second order model, incorporating finite Reynolds numbereffects in homogeneous flows. These models will be used in the analysis of our results. B. Experimental Method
In our investigation, we used the results from a wind-tunnel experiment. Detailed descriptions of theexperimental apparatus, measurements, and post processing, were given in Refs. [49, 50], yet for completeness,we present a brief overview here as well. We used the environmental wind-tunnel laboratory at IIBR, featuringa 14 meters long open wind-tunnel with a 2 × cross sectional area, that is compatible for conductingexperiments mimicking turbulent flows in the atmospheric surface layer. The canopy flow was modeled byplacing flat rectangular plates on the bottom floor of the wind tunnel. Our mixed height canopy layer wasconstructed of two types of plates with a height of either H or H , and a width of H , where H = 100 mm.The two types of plates were positioned in consecutive rows and at a staggered orientation, see Fig 1(a) and(b). The entire upstream part of the test section was fitted with roughness elements. The canopy frontal areaindex, defined as the plate frontal area divided by the lot area, Λ f = A f /A T = / , ( A f being is the elementfrontal area, and A T the lot area of the canopy). These densities categorizes our canopy as a moderatelydense and deep canopy.We gathered data at two levels of the free stream velocities, corresponding to Reynolds numbers Re ∞ ≡ U ∞ H/ν = 16 × and 26 × ; with U ∞ = 2 . − being the free stream mean velocity measuredwith a sonic anemometer at the center of the wind-tunnel cross section, and ν the kinematic viscosity of theair. In what follows we adopt the frame of reference commonly used in the canopy flow literature – x isstreamwise aligned longitudinally within the wind-tunnel, y is in the horizontal cross-stream direction, andthe positive z axis is directed vertically away from the bottom wall at which z = 0.In the experiment, we tracked fluid tracers using a real-time extension of the 3D-PTV method [49]. Theflow was seeded with hollow glass spheres of diameter 5 micrometer on average, with the Stokes number St = τ p /τ η ≈ .
05. The tracers were illuminated with a 10W, 532nm continuous wave laser; 3D positionswere inferred from 4 Mega pixel images at a resolution of 50 µ m per pixel and rates of 500 Hz within thecanopy layer z ≤ H , and 1000 Hz above the elements z > H . We implemented camera calibration, stereomatching and tracking [60], to reconstruct the tracer particle’s trajectories by using the OpenPTV [61] opensource software, integrated to operate with the real-time image analysis extension. The trajectory dataanalysis was performed by employing our open source Flowtracks package [62]. C. The Sub-Volume Approach
We analyzed trajectories in 20 sub-volumes centered at different locations. The sub-volumes are rectangularcuboids found at 4 horizontal locations, (see diagram in Fig. 1(a)), and at 5 different heights above the wind-tunnel bottom wall. The 4 horizontal positions of the sub-volumes are labeled alphabetically a , b , c and d ; a is found immediately downstream of a tall element, b is upstream to the next tall element, and c and d are parallel to the former elements and positioned around a short element. At each horizontal position, a − d , we used 5 vertical slabs of thickness δz = 0 . H which defines a total of 20 sub-volumes; the verticalslab position is labeled numerically 1-5 which correspond to the heights 0.5-0.7, 0.7-0.9, 0.9-1.1, 1.1-1.3 and1.3-1.5 H , respectively. Thus, for example, the sub-volume b . H < z ≤ . H . An animation of a sub-sample from our data set can be seen online through the link [63].Lagrangian statistics are commonly represented based on a common point of origin ( x , t ) and as a functionof time [54]. Thus the Lagrangian mean of an arbitrary function (cid:104) A ( x , t , t ) (cid:105) , is normally defined over flowensembles. In this work, we assume that the flow is stationary (see [49]) and ergodic, and thus we replace theabove with an average over t , presented as a function of time lag τ = t − t , where t is the time a particlewas first spotted in our measurement volume. In addition, we present Lagrangian statistics for particles with x in each of the sub-volumes, namely statistics are sub-volume averaged. Therefore, the ensemble averageover sub-volume V , (cid:104) A (cid:105) v of any quantity A , is calculated as: (cid:104) A ( x , τ ) (cid:105) v ≡ N N (cid:88) i =0 A i ( x v , τ ) (10)such that x v are positions of Lagrangian tracers inside the sub-volume x v ∈ V , and N is the number ofLagrangian trajectories in the sub-volume over the measurement period of time.It is notable that the sub-volumes we used are small. In particular, δz < (cid:104) u x (cid:105) ∂ (cid:104) u x (cid:105) /∂z , and we observed thatthe flow’s statistics had only minor variations within each sub-volumes. This frames the focus of our workon small scale Lagranigan statistics. D. Estimating the Mean Rate of Dissipation
The Eulerian second order structure function is defined as the second moment of spatial velocity differences δ u (cid:48) r ≡ u (cid:48) ( x ) − u (cid:48) ( x + r ) [54]. Assuming local isotropy and homogeneity, the Kolmogorov universal similaritytheory [64] predicts that the longitudinal component of the second order structure function, S LL , and itstransverse component, S NN , (i.e. the components aligned with r or normal to r respectively) admit to thefollowing scaling law in the inertial range of scales: S LL ( x , r ) = (cid:28)(cid:16) δ u (cid:48) r · r r (cid:17) (cid:29) = C ( (cid:15)r ) / S NN ( x , r ) = 12 (cid:16) (cid:104) δ v ( r ) (cid:105) − S LL ( r ) (cid:17) = 43 C ( (cid:15)r ) / (11)where in the homogeneous turbulence case the x dependence drops, r = | r | , and C ≈ . S LL in each sub-volume using the Lagrangiandataset and averaging velocity differences over spherical shells. This gives a structure function that is isotropicby construction, namely it depends only on the distance r . For example, in the main panel of Fig. 1(c) wepresent our estimations of S LL ( r ) in the sub-volumes b1, b3 and b5 for the Re ∞ = 16 × case. Note thatour estimation of S LL does not use the Taylor’s hypothesis.Using S LL we estimated an empirical mean rate of dissipation in each sub-volume. Since S LL is quadraticwith r in the dissipation range and should change very slowly at large scales above the integral scale L [54],FIG. 1: (a) A top view sketch of the canopy repeating unit cell. Grey shaded regions show the positions ofthe four sub-volumes. Black thin rectangles represent the canopy roughness elements. (b) An isometricsketch of the short and tall roughness obstacles used. (c) Main panel shows the second order Eulerianlongitudinal structure function in three sub-volumes, and dashed lines show the isotropic model, Eq. (11),with the estimated values (cid:15) v . The inset shows the compensated structure functions for the same cases. (d)The ratio between the transverse and the longitudinal Eulerian structure functions, compared with the K41isotropic value 4 / S LL ( r ) /C ( (cid:15) r ) / should peak at an intermediate range η < r < L , andthis is shown in the inset of Fig. 1(c). Thus, we defined a sub-volume averaged dissipation rate as (cid:15) v ≡ max r (cid:20) S LL ( r ) C r / (cid:21) (12)where S LL is estimated with samples x ∈ V . Using Eq. (12) instead of a least square fitting has the advantageof not having to specify a range of r where an inertial range scaling supposedly exists; on the other hand, itmay overestimate the value that would have been obtained in a fitting process, but this uncertainty is low(arguably much lower than the uncertainty in the value of C ). Also, due to anisotropy and inhomogeneity,the applicability of Eq. (11) to canopy flow turbulence is not straightforward, as observed for example inRefs. [23, 66, 67]. Nevertheless, in Fig. 1(c) we compare our estimations of S LL ( r ) with Eq. (11) using (cid:15) v ,and for the three cases we observe a distinct range of r , in which an agreement between the theory andthe empirical data exists. In addition to that, the transverse components, S NN were similarly estimatedusing Eq. (11) and the Lagrangian dataset. The ratios S NN /S LL were calculated and an example is shownin Fig. 1(d) as a function of r for the three sub-volumes b1, b3, and b5. The K41 isotropic value of S NN = S LL [64] is plotted for comparison as well. The figure shows that for all three sub-volumes, theratio S NN /S LL decreases with r while crossing the 4/3 value at a limited range of r/η . These observationsof agreement with K41 similarity in a certain limited range of r supports the use of (cid:15) v as a parameterizationin our experiment. Therefore, sub-volume averaged dissipation scales were calculated as η ≡ (cid:0) ν /(cid:15) v (cid:1) / and τ η ≡ ( ν/(cid:15) v ) / being the length and time scales respectively, in analogy to the usual case [68].With the estimations of (cid:15) v , a Taylor microscale Reynolds number is defined as – Re λ,v ≡ v (cid:48) (cid:18) (cid:15) v ν (cid:19) / (13)where v (cid:48) = (cid:104) (cid:80) i ( v i − (cid:104) u i (cid:105) ) (cid:105) is the RMS of particle velocity relative to the sub-volume averaged velocity.Estimations of (cid:15) v and Re λ,v were repeated for all sub-volumes and the values are tabulated in the Appendix.The Reynolds numbers Re λ,v varied with x in the range 350-850 due to the flow’s inhomogeneity, and hadnegligible dependence on Re ∞ . The obtained values were used for parameterization of the results presentedbelow and the subscript “ v ” is omitted to facilitate readability. III. QUASI-HOMOGENEOUS LAGRANGIAN DYNAMICS IN THE SMALL SCALES
In this section, we determine the dominant dynamical factors in the context of our measurements byquantifying and comparing contributions from the different terms in the LSM eq. (7). Identifying the mostrelevant processes in the context our Lagrangian measurements is crucial since it will present our results inthe proper context. In addition to that, it will dictate the tools that we use in the consecutive analysis.For the purpose of this analysis, we use the Thomson formulation [19] LSM for inhomogeneous turbulence,eq. (1) and eq. (7). Although the LSM is an equation of motion in the Lagrangian framework, its inputis a set of Eulerian variables, so to use it we must cast our Lagrangian measurements onto an Euleriancoordinate system. Specifically, we require field representations of the mean velocity, (cid:104) u i (cid:105) , the turbulent stresstensor, R ij , and the Lagrangian structure function parameterization, C (cid:15) , all of which are available from ourempirical dataset. Thus, we estimate the velocity statistics by using the sub-volume averages: (cid:104) u i (cid:105) = (cid:104) v i (cid:105) and R ij = (cid:104) v (cid:48) i v (cid:48) j (cid:105) . For C (cid:15) , we rely on the fact that through our Lagrangian measurements we can calculate D ij ( τ ) directly from the definition, Eq. (5). This is unlike Eulerian measurements that must rely on modelsor on previous measurements in similar flows. Thus, we calculate the structure function parameterizationby C (cid:15) = max τ [ D ii τ ]; this issue is discussed in detail below. Following that, we obtain continuous anddifferentiable field estimations by using inverse-distance weighted interpolations of the sub-volume averageddata (for details, see Appendix D). Through this interpolation scheme, we obtained estimations of the Eulerianfields in our volume of measurement, that allow to calculate the models’ coefficients.We compare the magnitudes of the various LSM terms in Eq. (7) along the trajectories from our experimentusing the above interpolation schemes . In particular, to facilitate the analysis, we employ the Wilson etal. [56] formulation for the models coefficients: T i = 12 ∂R ij ∂x j ; T ijk = 12 R − lj ∂R il ∂x k ; T ij = T ijk (cid:104) u k (cid:105) In the top panel of Fig. 2 we show probability distributions for the magnitude of the various LSM termsalong the trajectories in our empirical dataset in sub-volume b3, namely at height 0 . < zH ≤ .
1. The curvesessentially show the magnitudes of the forces that acted on our measured particles, decomposed according tothe LSM. The T i term accounts for inter sub-volume variations of the Reynolds stress tensor, and it has thenarrowest distribution of all the terms. The T ij and T ijk terms are characterized by much wider PDFs sincethey include both spatial variations of R ij and the randomness of temporal velocity fluctuations. The figurealso shows that the term accounting for relaxation of turbulent fluctuations due to dissipation, C (cid:15)R − ij , isthe largest of all the terms. This implies that velocity changes due to effects of random turbulent fluctuationswere generally more dominant than the effects of all the other force components. Notably, the T m termsarise due to inhomogeneity of the flow, while the C (cid:15)R − ij term accounts for turbulent fluctuations thatexsist in both homogeneous and inhomogeneous flows (Eq. (2)). Thus, in the bottom panel of Fig. 2 wecompare PDFs of | φ/g | = | T i + T ij v (cid:48) j + T ijk v (cid:48) j v (cid:48) k | that represents the magnitued of the combined effect offlow inhomogeneity, with that of the dissipation term. The data is shown for three sub-volumes representingthe regions inside, right at top of, and above, the canopy layer (b1, b3, and b5 respectively). The mean ofeach PDF is also marked with a vertical line. The figure shows that for all the regions tested the relaxationdue to random fluctuations were more dominant than the effect of flow inhomogeneity. In particular, theratio between the mean value of the terms was (cid:104) | C (cid:15)R − ij v (cid:48) j | | φ i /g | (cid:105) ≈ . | φ i /g | was generally smaller than the homogeneous term, | C (cid:15)R − ij v (cid:48) j | . This means that randomforces due turbulent dissipation were stronger than forces due to flow inhomogeneity. In regions fartherabove the canopy layer (i.e. sub-volumes 4 and 5 with z > . H ), the difference between the two termsFIG. 2: Probability distribution functions (PDFs), where the horizontal axis stands for magnitude ofvarious inhomogeneous LSM terms shown in the legend, sampled along the Lagrangian trajectories from ourdataset. Data corresponding to all the model’s terms in sub-volume b3, representing 0 . H < z ≤ . H , isshown in the top panel. PDFs for the sum of flow inhomogeneity contributions ( φ i /g ) in three sub-volumes,representing regions inside (b1), at the top of (b3), and above (b5), the canopy layer, is shown in the bottompanel. Vertical dotted and dashed lines mark the mean values of the PDFs. For all cases Re ∞ = 1 . × .reduced, although the turbulent fluctuations through | C (cid:15)R − ij v (cid:48) j | were still significantly more dominant.In comparing the variation of the terms with z we observed that | C (cid:15)R − ij v (cid:48) j | typically increased with z inside the canopy, peaked at the top of the canopy, and than reduced with z above the canopy; on the otherhand, both | T ij v (cid:48) j | and | T ij v (cid:48) j | became increasingly stronger with the height. Thus, as the distance from thewall increases the effects of dissipation became weaker and those of inhomogeneity became stronger. Thissuggests that the reason for this effect is the direct interaction of the flow with canopy obstacles, that itseffects are weaker above the canopy. The flow-obstacle interaction is known to cause an increased energyproduction and dissipation due canopy drag and wake production [17], and seems to dominate the effects ofspatial variations of flow statistics in the canopy obstacles’ wake region.As described above, our analysis of Lagrangian statistics was performed in sub-volumes of space that weremuch smaller than the integral scale of the flow (i.e. of the order H ), and indeed the flow statistics didnot change appreciably inside each sub-volume. Furthermore, the results presented in this section showunequivocally that the flow inhomogeneity did not have a dominant effect on the Lagrangian dynamics. Thisleads to the conclusion that there is a quasi-homogeneous regime of Lagrangian statistics at small times, sincethe contributions from flow inhomogeneity are negligible at small scales as compared to turbulent dissipation.It is important to stress that the canopy flow is inherently inhomogeneous since, for example, statistics changefrom one sub-volume to another. Nevertheless, inhomogeneity effects on short time Lagrangian statistics, i.e.in scales (cid:28) (cid:104) u x (cid:105) ∂ (cid:104) u x (cid:105) /∂z , were sub-dominant as compared to the strong effects of dissipation. Our observationthus reveal the existence of a quasi-homogeneous regime at short-times and small-scales. IV. DIRECT ESTIMATION OF LAGRANGIAN STATISTICS
In this section, we use the wind-tunnel Lagrangian dataset in order to extract two critical LSM parameters:the Lagrangian velocity decorrelation timescale and the Kolmogorov constant, C . To do so we calculate theautocorrelation and the structure functions and directly extract the two parameters from their definitions. A. Lagrangian Autocorrelation and Decorrelation Timescale
We estimated the Lagrangian autocorrelation functions using a formula equivalent to Eq.(2.6) in the paperby Guala et al. [37] with the ensemble averaging Eq. (10), as presented and discussed in Appendix B. Fig. 3(a)presents ρ ii for the three components of Lagrangian velocity in sub-volume b3. For all three components, aconcave shape is seen at the origin, reminiscent of a parabolic decrease at small delay times, providing theway to estimate the Taylor micro-timescale [69]. The autocorrelation functions decrease monotonically withthe increasing time lag, τ . The rate of decrease is roughly the same for ρ xx and ρ yy , whereas the decreaseis faster for the vertical component, ρ zz . These observations were robust throughout all of the sub-volumesand the two Re ∞ cases.The Lagrangian autocorrelation function ρ ii ( τ ) does not decrease to zero within the range of our mea-surements, so we cannot use the integral in Eq. (3) to estimate the T L,i directly. Therefore and similarly toprevious Lagrangian measurements (Refs. [34, 70]), we define a Lagrangian decorrelation timescale , T i , thatwe obtain by fitting our results to the autocorrelation function of an LSM. In accordance with the results ofSection III, we use here a LSM for homogeneous flows. Specifically, we used a least square minimization tofit our measurements to the autocorrelation function in Sawfords model, eq. (9), to obtain T ,i and T ,i ineach sub-volume and then define T i = T ,i + T ,i . It is noted that a calculation of T L,i using the full integralof ρ ii according to the theory, may result in a larger timescale on the order of the turnover timescale for largecoherent structures above the canopy. The concave shape of the autocorrelation function at τ → T i for the three velocity components by fitting the data as shown in Fig. 3(a). The fit rangewas limited to the time lags that correspond to a half of the sub-volume crossing time, in order to avoidthe possible finite volume effects (as discussed in Appendix C and Ref. [71]). This is a common approach inexperimental data analysis because of increased uncertainty of correlations at larger time lags [72]. Lastly,we note that every data point of ρ ii corresponds to the average of at least 15 × samples, where the relativemean squared error was of the order of a few percents.It is worth noting that the Lagrangian velocity become decorrelated, namely ρ ii ( τ ) reduced considerably,while the particles were still within the small sub-volumes we used. This is in agreement with an observationin our previous paper [49] where we detected the Taylor asymptotic dispersion regime [73], and also withour observation in Section III of the quasi-homogeneous regime.The empirical Lagrangian decorrelation times, T i , are presented in Fig. 3(b,c,d) for each velocity componentas a function of height for the two Reynolds numbers tested. Inside the canopy layer, the values of T x areroughly constant and they increase above the canopy. In contrast, T y values are highest at the lowest sub-volume yet retain a roughly constant value above the canopy. The values of T z are the lowest of the threecomponents and show only minor variation with height above H . Furthermore, the Lagrangian decorrelationtimes are consistently higher for the Re ∞ = 16 × case, as compared to the higher Re ∞ = 26 × case.The distributions of T i can be associated with physical processes that are known to occur in canopy flows.The increase of T x above the canopy may be related to large scale coherent structures that are known toexist above canopies due to the shear instability, aka the mixing layer analogy [17, 24]. The increase of T y inside the canopy layer is attributed to the change in the roughness density with height – the lower roughnesselements caused an increased frontal area density and increased the ”shielding” (e.g. [74]), which contributedto a tunneling effect of a cross-flow inside the canopy layer. Lastly, lower values of T z as compared to T x and T y are in agreement with estimations of the Eulerian integral timescales from velocity measurement, forinstance by Refs. [18, 43, 75], which may be associated with an inclined orientation of coherent structuresthat was reported in the literature, e.g. by Shaw et al [75] using two-point Eulerian correlations.0 .
00 0 .
02 0 .
04 0 . τ [s]0 . . . . . . ρ ii ( τ ) Fit Range xyz (a) .
000 0 .
025 0 .
050 0 .
075 0 .
100 0 . T x [s]0 . . . . . z / H (b) .
000 0 .
025 0 .
050 0 .
075 0 .
100 0 . T y [s]0 . . . . . z / H (c) (d) FIG. 3: (a) Autocorrelation function of Lagrangian velocity in the sub-volume b Re ∞ = 16 × ,presented against the time lag, shapes denote data points and lines represent fits to the model Eq. (9).(b)-(d) The Lagrangian timescale for the x , y and z velocity component as a function of sub-volume height.Filled shapes stand for Re ∞ = 16 × and empty shapes for Re ∞ = 26 × . B. Second Order Lagrangian Structure Function
Using the sub-volume averaging, we estimated the Lagrangian second order structure function through itsdefinition, Eq. (5). The results for trajectories in sub-volume b3 are shown in Fig. 4(a) in a compensatedform, D ii /(cid:15)τ , which corresponds to the Kolmogorov scaling in the inertial range. For high Reynolds numberHIT flows, the existence of Kolmogorov scaling in a Lagrangian inertial range would lead to a plateau inthe compensated plot. The figure shows that such a plateau does not appear in our data, and instead, onlynarrow peaks are seen. Such peaks are characteristic of low to moderate Reynolds number flows, and similarobservations were reported in numerous previous studies involving other types of flows, Refs. [33, 36, 40, 70,71, 76–79]. Therefore, we use the typical empirical estimate of C ,i , which is defined using the height of thepeaks C ,i = max τ (cid:20) D ii ( τ ) (cid:15) τ (cid:21) . (14)1 τ /τ η D ii ( τ ) (cid:15) τ D xx D yy D zz (a) C . . . . . z / H abcd (b) FIG. 4: (a) Second order Lagrangian structure function of velocity differences in sub-volume b Re = 16 × , presented as a function of time normalized with the inertial-range scaling. The continuous,dashed, and dash-dotted lines represent Sawford’s second order LSM [58] fitted to the x , y and z components respectively. (b) The estimated Kolmogorov constant Eq. (14) from the x velocity component,for all sub-volumes and two Reynolds numbers as a function of height. Full shapes for Re ∞ = 16 × andhollow for Re ∞ = 26 × .In the case of sub-volume b3, shown in Fig. 4(a), the values obtained are in the range C ,i ∈ (6 . , .
0) for thethree velocity components, meaning weak anisotropy ∼ O (10%).We applied Eq. (14) to the D xx ( τ ) obtained in the different sub-volumes in the canopy flow model. Theobtained values of C ,x are presented in Fig. 4(b) as a function of height. Values of C ,i calculated with D yy and D zz were very similar to Fig. 4(b) and are not shown for the sake of brevity. A considerable scatter isseen in the values, roughly ∼ Re λ variations, or due to a sensitivity of Eq. (14) to small uncertainties in the structure-function.In Fig.4(b), C ,x does not show a dependence on Re ∞ , which is consistent with the behavior of Re λ . Goingfrom inside the canopy and increasing in z , C ,x initially increases and roughly levels off at C ,x ≈ . . < z ≤ .
2. Further up above the canopy, C ,x decreased even though Re λ was seen to be highestat this height level.The non-monotonous behavior of C ,x with z suggests an important conclusion. In HIT, C is governed by Re λ alone (e.g. as suggested in [78, 80]), however, the fact that here C is not monotonous with z while Re λ increases with z monotonously suggests that in the canopy flow C depends on other parameters in additionto Re λ . An explanation for this observation can be offered through an analysis by Poggi et al. [18]: theeffect of wake production on the Lagrangian structure function, leading to scale dependence of the rate ofdissipation in canopy flows (i.e. the spectral bump [17]), is ”lumped” into C . Therefore, according to thisexplanation, the observed increase of C at the top of the canopy in our measurements is due to a strongwake production, injecting turbulent kinetic energy at small scales.In Section IV A it was observed that the autocorrelation function was concave at the origin (Fig. 3), andthat we obtained a good fit for the data using the second-order LSM, Eq. (8). In a homogeneous flow, theLagrangian structure function and the D ii are simply related by D ii ( τ ) = 2 σ i [1 − ρ ii ( τ )] , (15)and thus with regard to the observed quasi-homogeneity of our flow it would be instructive to examinethis relation here as well. Thus, we used Eq. (15) to fit the empirical data for D ii ( τ ), where we used theexpression for ρ ii ( τ ) from the second order LSM, Eq. (9). Specifically, we used a least square algorithm tofit the three model parameters for the three components of the structure function. The resulting curves areshown with lines in Fig. 4(a). The good match that was obtained for the empirical data shows that single-particle statistics in our measurements are represented well by Sawford’s second order LSM [58], Eq. (8).2Eq.(17)(b)(a)FIG. 5: Histograms of Lagrangian de-correlation times in the x direction from all sub-volumes and two Re ∞ (a total of 40 points), normalized with the (a) Eq. (16), and (b) the dissipation timescale.This further reinforces the picture of quasi-homogeneity of Lagrangian statistics at short times in our flow. V. LAGRANGIAN RAPID DECORRELATIONA. Observation of Rapid Decorrelation in the Canopy Flow
In light of the small-scales’ quasi-homogeneity, we find it instructive to compare our empirical estimates of T i with previous results for purely homogeneous turbulent flows. We first take the LSM Eq. (5), accordingto which T L,i = 2 v (cid:48) i C ,i (cid:15) . (16)This relation was previously used to estimate the Lagrangian integral timescale based on Eulerian measure-ments in canopy flows, for example in Refs. [12, 15, 23]. The comparison of our T i with Eq. (16) is shownin Fig. 5(a) through a histogram of the property T x / ( σ u C (cid:15) ), taking values from all the sub-volumes and two Re ∞ . The histogram shows a large scatter of values in the range [0.3–1] with an average of 0.65, and thusimplies that the empirically estimated T i is significantly shorter than Eq. (16).The second comparison, in Fig. 5(b), compares the separation of scales T i /τ η in our canopy flow with thatof HIT at similar Re λ . A histogram of the decorrelation times, T x , normalized with the dissipation timescale, τ η , is shown in Fig. 5(b). The results fall into two groups with values in the order of T x ∼ − τ η , and T x ∼ − τ η that were seen to occur at different height levels, consistent with the increase of Re λ fartheraway from the wall. The values of T x /τ η seen in Fig. 5(b) are low as compared to values of T L /τ η usuallyencountered in HIT at comparable Re λ . For example, Sawford et al. [81] suggested T L τ η = (cid:34) .
77 + (cid:18) Re λ . (cid:19) / (cid:35) / (17)following the empirical fit to DNS data. Using the Re λ values from our measurements, the estimates basedon Eq. (17) are plotted in the inset of Fig. 5(b). Here the values are seen to be roughly an order of magnitudehigher than those obtained by directly fitting the autocorrelation functions. Therefore, Fig. 5 shows that theseparation of scales in our canopy flow is much smaller than what would have been expected in comparableHIT case.3In Fig. 5, T i was compared with properties estimated using spatial information of the flow, namely (cid:15) . Wecan reinforce the above observation of relatively small separation of scales by using a purely Lagrangianproperty that characterizes the small scales, such as the particles’ acceleration, a = d v dt . Thus, we contrastthe autocorrelation of v i and a i in our canopy flow with those in a comparable HIT flow. For the HITflow, we use the DNS data available from the Johns Hopkins Turbulence Database (JHTDB) [79, 82]. Wedownloaded this benchmark dataset in a previous study [46], and we use it here again to compare theautocorrelation functions. Autocorrelation functions from the DNS data and the data from sub-volume b3are plotted together in Fig. 6. Note that the Re λ in both cases are very similar, Re λ ≈
440 in sub-volume b3and Re λ ≈
433 in the DNS [79], and so they should present similar separation of scales if T L /τ η = f ( Re λ ).Fig. 6 shows that the autocorrelation of the acceleration components decay at rates very similar to each otherin both cases. However, the velocity in the canopy flow decorrelates much faster than in the HIT case. Tobe more quantitative, denoting by T a the time of the first zero crossing of the acceleration, it was found that T x /T a ≈ . T x /T a ≈ . T a was roughly the same in both cases, the velocity became decorrelated roughly 3 times faster in thecanopy flow as compared to the HIT case. This observation reinforces the observation that the separation oftimescales in our canopy flow is smaller than in the comparable HIT case. b3b3 Ta T x ≈ . T a T x ≈ . T a FIG. 6: Lagrangian auto correlation functions of the acceleration component and the velocity component ofLagrangian particles. Hollow symbols correspond to trajectories for the DNS [46, 79, 82] and filled symbolsto canopy trajectories from sub-volume b3 at Re ∞ = 16 × .Together, Figs. 5 and 6 demonstrate that the Lagrangian velocity components in our canopy flow becamedecorrelated much faster than what would have been expected in a homogeneous isotropic turbulent flows ata similar Reynolds number. In particular, the above demonstrates that in canopy flows, unlike the HIT case, T i /τ η is not a function of the Reynolds number alone. This observation will be termed in what follows rapiddecorrelation . B. Turbulence–Obstacle Interaction as the Source for Rapid Decorrelation
In in section III, we observed a dominance of turbulent fluctuations over inhomogeneity in the Lagrangiandynamics, and in section IV B we suggested that the structure function constant C was affected by the wakeproduction in accordance with the arguments of Poggi et al. [18]. These observations show a strong influenceof small scale dynamics on Lagrangian statistics in our measurements. These considerations lead us to putforth the notion that the observed rapid-decorrelation was also a consequence of wake effects due to the directobstacle-flow interaction. In this section we examine this conjecture.The encounter of the the flow with canopy obstacles leads to generation of drag that injects turbulent kineticenergy at flow scales with sizes that are determined by the geometry of the roughness obstacles, so-calledwake-production [17, 83]. Therefore, if our conjecture was true, we would expect that the dispersion in thewakes will be dominated by flow disturbances with a similar size. The appropriate length scale of dispersion4 (a) (b) FIG. 7: (a) Lagrangian dispersion length scale L x = T x σ x , presented as a function of height, both axesnormalized by H . (b) Two-particle spatial correlation function, plotted against distance normalized by thecanopy top height. Data for particles in sub-volume b3.is L i = σ i T L,i , which was observed to be correlated with the scale of forcing in previous experiments(Refs. [70, 84]). The width of our obstacles was d = 0 . H , so if the hypothesis is true we should see L i ∼ d .Thus, we calculated L x = T x σ x in the different sub-volumes and the results are presented in Fig. 7(a) againstheight. Inside the canopy layer, z ≤ H , L x is nearly constant, L x ≈ . H = 2 . d independently of Re ∞ ,which may well be due to such wake disturbances. Above the layer, L x increases, reaching ≈ . H at thehighest sub-volume. The lack of deviation of L x inside the canopy along with increase of L x at z > H areconsistent with the notion of weakening of the wake’s influence above the canopy. Therefore, the estimatedvalues of L x are consistent with the conjecture that rapid decorrelation occurred due to the flow disturbancesin the obstacles’ wakes.To further support our conjecture, we wish to demonstrate that a disturbances occurred in our flow at thescale L x . We demonstrate this through the two-particle spatial velocity correlation function. Specifically, letus define ρ (1 , i ( r ) = (cid:104) v (cid:48) (1) i ( t ) v (cid:48) (2) i ( t ) | r ( t ) (cid:105) Var[ v (cid:48) i ( t )] r (18)where v (cid:48) (1) i and v (cid:48) (2) i are the velocity fluctuation of two different particles at the same time, and where theaverage in the numerator is performed for pairs of trajectories with a distance of r ( t ) = | x (1) ( t ) − x (2) ( t ) | between them. In the denominator, Var[ v (cid:48) i ( t )] r is the variance of the velocity components calculated overthe same ensemble of particles used in the numerator. Note that ρ (1 , i is a correlation with no separation intime but only in space, and thus it can be used to examine the spatial structure of the flow. Also note that ρ (1 , i ( r ) is analogous to the Eulerian two-point spatial velocity correlation function (for example see Shaw etal. [75]), however, ρ (1 , i ( r ) is isotropic by construction since the average is performed over spherical shells.The two-particle correlation of v (cid:48) x , calculated using trajectories from sub-volume b3, is shown in Fig. 7(b).The same data is shown in linear-log scales in the main figure and in linear scales in the inset. The ρ (1 , i ( r )decreases monotonously with r . In the range r (cid:46) . H , the correlation decreases faster than at r (cid:38) . H .In the linear-log scales, the data at each interval roughly fits a straight line, where in each interval it has adifferent slope. Therefore, the data points were fitted with exponential decays with a different rate in each ofthese two ranges of r , which provided a good approximations of the data. The two fits, shown in dashed anddot-dashed lines, highlight the transition of the ρ (1 , i ( r ) from one rate of decay to another that occurs rightat r ≈ . H ≈ L x . This transition of ρ (1 , i ( r ) from one rate of decay to another at r ≈ L x may suggest anexistence of flow disturbances of characteristic size L x ∼ d , namely corresponding to the width of the flowobstacles. Furthermore, such a transition of ρ (1 , i ( r ) at r ≈ L x was robust for the sub-volumes inside the5canopy, however not above it. Thus, since T x = L x /σ x , we can speculate that such disturbances at the wakescale may have lead to the observed rapid-decorrelation; nevertheless we cannot prove this conjecture at thistime.The two pieces of evidence presented above are in agreement with our conjecture, and thus leave the notionthat wake production is the main cause of the rapid decorrelation of Lagrangian velocity a valid possibility. Aconclusive proof will require further exploration, for example, by using a flow with various degrees of spectralshort-circuiting versus inhomogeneity effects. VI. DISCUSSION & CONCLUSIONS
In this work, we used experimental measurements to estimate Lagrangian statistics in a wind-tunnel canopyflow model directly in the Lagrangian framework. Our analysis indicates that turbulence-obstacle interaction,through wake-production had a significant effect on the short time Lagrangian statistics in our flow, and isrelevant for Lagrangian stochastic models. In particular, our key result is that in spite of the large scaleinhomogeneity (e.g. Figs 3, 4(b), 5 and 7(a)), we detected a quasi-homogeneous regime of Lagrangianstatistics at short times and on the small scales. Furthermore, we show that the spatial variations of theseparation of scales, T i /τ η , and the Kolmogorov constant, C , cannot be explained by the variation of theReynolds number, Re λ ; this suggests that unlike in HIT, they depend on additional parameters other than Re λ . The main difference is that the decorrelation timescale of the Lagrangian velocity is much shorter thanwould have been expected in a homogeneous case. We thus infer that both the rapid decorrelation and thealteration of C are direct consequences of strong wake production.The strong influence of the wakes on Lagrangian dynamics had important implications on our analy-sis. First, we found that the Lagrangian statistics in the quasi-homogeneous regime are recovered wellby the second order LSM for homogeneous flows. Second, due to the small separation of scales, so-calledrapid-decorrelation, we detected significant finite Reynolds number effects on the Lagrangian autocorrelationfunctions, and this is despite the fact that Re λ in our canopy flow was rather high, in the range of 350–850.Essentially, this is a demonstration that finite Reynolds number effects can be important in cases where theRichardson-Kolmogorov cascade is short-circuited.We expect that our results will be relevant for modeling short range dispersion in flows with intensedissipation and spectral short-circuiting (in particular, where | C (cid:15)R − ij v (cid:48) j | (cid:29) | φ i /g | ). We achieved thisthrough high frontal area density (i.e. λ f = / ) leading to strong drag, thin obstacles that producedturbulent kinetic energy at a rather small scale, and obstacles with variable heights in consecutive rows.To conclude, the observations presented in this work make up a unique view on the Lagrangian dynamicsin the canopy flows in the small scales. It is our view that short term dispersion modeling in canopy flowsthrough LSMs may achieve increased accuracy by paying particular attention to wake dynamics in the canopyflows. Our work also highlights the importance of gathering Lagrangian statistics directly in the Lagrangianframework that is becoming possible with recently introduced technologies (i.e. [49] and references therein).Other important topics that were not dealt with here include the effects of mechanical diffusion, the mixing-layer analogy, and the inhomogeneity, on Lagrangian statistics, that due to the small scale of our observationcould not have been assessed here, and thus leave considerable scope for future investigations. ACKNOWLEDGEMENT
We are grateful to Meni Konn, Sabrina Kalenko, Gregory Gulitski, Valery Babin, Amos Shick andMordechai Hotovely, for their assistance in preparing and performing the wind tunnel experiment. Wealso thank an anonymous reviewer for provoking an important discussion. This study is supported by thePAZY grant number 2403170. [1] R. E. Britter and S. R. Hanna. Flow and dispersion in urban areas.
Annual Reviews in Fluid Mechanics ,35:469–496, 2003. [2] R. Nathan, G. G. Katul, H. S. Horn, S. M. Thomas, R. Oren, R. Avissar, S. W. Pacala, and S. A. Levin.Mechanisms of long-distance dispersal of seeds by wind. Nature , 418, 2002.[3] M. R. Raupach and A. S. Thom. Turbulence in and above plant canopies.
Annual Review of Fluid Mechanics ,13(1):97–129, 1981.[4] M. R. Raupach. A lagrangian analysis of scalar transfer in vegetation canopies.
Quarterly Journal of the RoyalMeteorological Society , 113(107-120), 1987.[5] E. Gavze and E. Fattal. Description of a turbulent cascade by a Fokker-Planck equation.
Boundary LayerMeteorology , 169:297–326, 2018.[6] M. R. Raupach. Applying lagrangian fluid mechanics to infer scalar source distributions from concentrationprofiles in plant canopies.
Agricultural and Forest Meteorology , 47:85–108, 1989.[7] T. K. Flesch and J. D. Wilson. A two-dimensional trajectory-simulation model for non-gaussian, inhomogeneousturbulence within plant canopies.
Boundary Layer Meteorology , 61:349–374, 1992.[8] J.D. Wilson and B.L. Sawford. Review of lagrangian stochastic models for trajectories in the turbulent atmo-sphere.
Boundary-Layer Meteorology , 78:191–210, 1996.[9] M. W. Rotach, S.-E. Gryning, and C. Tassone. A two-dimensional lagrangian stochastic dispersion model fordaytime conditions.
Quarterly Journal of the Royal Meteorological Society , 122(530):367–389, 1996.[10] D. Baldocchi. Flux footprints within and over forest canopies.
Boundary-Layer Meteorology , 85(2):273–292, Nov1997.[11] R. Leuning, O. T. Denmead, A. Miyata, and J. Kim. Source/sink distributions of heat, water vapour, carbondioxide and methane in a rice canopy estimated using lagrangian dispersion analysis.
Agricultural and ForestMeteorology , 104:233–249, 2000.[12] D. E. Aylor and T. K. Flesch. Estimating spore release rates using a lagrangian stochastic simulation model.
Journal of Applied Meteorology , 40(7):1196–1208, 2001.[13] R. W. Arritt, C. A. Clark, A. S. Goggi, H. Lopez Sanchez, M. E. Westgate, and J. M. Riese. Lagrangian numericalsimulations of canopy air flow effects on maize pollen dispersal.
Field Crops Research , 102(2):151 – 162, 2007.[14] Gleicher S. C., Chamecki M., Isard S. A., Pan Y., and Katul G. G. Interpreting three-dimensional spore concen-tration measurements and escape fraction in a crop canopy using a coupled eulerian–lagrangian stochastic model.
Agricultural and Forest Meteorology , 194:118 – 131, 2014.[15] T. Duman, A. Trakhtenbrot, D. Poggi, M. Cassiani, and G. G. Katul. Dissipation intermittency increases long-distance dispersal of heavy particles in the canopy sublayer.
Boundary-Layer Meteorology , 159(1):41–68, Apr2016.[16] E. Fattal, O. Buchman, and E. Gavze. A lagrangian-stochastic model for pollutant dispersion in the urbanboundary layer over complex terrain - comparison with haifa campaigns. , 8A.4, 2018.[17] J. Finnigan. Turbulence in plant canopies.
Annual Review of Fluid Mechanics , 32:519–571, 2000.[18] D. Poggi, G. G. Katul, and M. Cassiani. On the anomalous bahavior of the lagrangian structure function similarityconstant inside dense canopies.
Atmospheric Environments , 2008.[19] D. J. Thomson. Criteria for the selection of stochastic models of particle trajectories in turbulent flows.
Journalof Fluid Mechanics , 180:529–556, 1987.[20] A. M. Reynolds. A second-order lagrangian stochastic model for particle trajectories in inhomogeneous turbulence.
Quarterly Journal of the Royal Meteorological Society , 125:557, 1999.[21] W. J. Massman and J. C. Weil. An analitical one-dimensional second-order closure model for turbulence statisticsand the lagrangian time scale within and above plant canopies of arbitrary structure.
Boundary-Layer Meteorology ,91:81–107, 1999.[22] ¨U Rannik, M. Aubinet, O. Kurbanmuradov, K. K. Sabelfeld, T. Markkanen, and T. Vesala. Footprint analysisfor measurements over a heterogeneous forest.
Boundary-Layer Meteorology , 97(1):137–166, Oct 2000.[23] D. Poggi, G. Katul, and J. Albertson. Scalar dispersion within a model canopy: Measurements and three-dimensional lagrangian models.
Advances in Water Resources , 29(2):326 – 335, 2006. Experimental Hydrology:A Bright Future.[24] M. R. Raupach, J. J. Finnigan, and Y. Brunet. Coherent eddies and turbulence in vegetative canopies: Themixing-layer analogy.
Boundary-Layer Meteorology , 78:351–382, 1996.[25] R. H. Shaw and I. Seginer. Calculation of velocity skewness in real and artificial plant canopies.
Boundary LayerMeteorology , 39:315–332, 1987.[26] H. M. Nepf. Drag, turbulence, and diffusion in flow through emergent vegetation.
Water Resources Research ,35(2), 1999.[27] S. B. Pope and Y. L. Chen. The velocitydissipation probability density function model for turbulent flows.
Physicsof Fluids A: Fluid Dynamics , 2(8):1437–1449, 1990.[28] G. A. Voth, K. Satyanarayan, and E. Bodenschatz. Lagrangian accelertion measuremetns at large reynoldsnumbers.
Physics of Fluids , 10(9):2268, 1998.[29] S. Ott and J. Mann. An experimental investigation of the relative diffusion of particle pairs in three-dimensionalturbulent flow.
Journal of Fluid Mechanics , 422:207–223, 2000. [30] N. Mordant, J. Delour, E. L´eveque, A. Arn´eodo, and J.-F. Pinton. Long time correlations in lagrangian dynamics:A key to intermittency in turbulence. Phys. Rev. Lett. , 89:254502, Dec 2002.[31] L. Biferale, G. Boffetta, A. Celani, B. J. Devenish, A. Lanotte, and F. Toschi. Lagrangian statistics of particlepairs in homogeneous isotropic turbulence.
Physics of Fluids , 17, 2005.[32] M. Bourgoin, N. T. Ouellette, H. Xu, J. Berg, and E. Bodenschatz. The role of pair dispersion in turbulent flow.
Science , 311:835–838, 2006.[33] P. K. Yeung, S. B. Pope, and B. L. Sawford. Reynolds number dependence of lagrangian statistics in largenumerical simulations of isotropic turbulence.
Journal of Turbulence , 7:N58, 2006.[34] N. Ouellette, H. Xu, M. Bourgoin, and E. Bodenschatz. Small-scale anisotropy in lagrangian turbulence.
NewJournal of Physics , 8, 2006.[35] J. Berg, B. Luthi, J. Mann, and S. Ott. Backwards and forwards relative dispersion in turbulent flow: anexperimental investigation.
Physical Review E , 74, 2006.[36] J. Bec, L. Biferale, M. Cencini, A. S. Lanotte, and F. Toschi. Effects of vortex filaments on the velocity of tracersand heavy particles in turbulence.
Physics of Fluids , 18(8):081702, 2006.[37] M. Guala, A. Liberzon, A. Tsinober, and W. Kinzelbach. An experimental investigation on lagrangian correlationsof small-scale turbulence at low Reynolds number.
Journal of Fluid Mechanics , 574:405–427, 2007.[38] R. J. E. Walpot, C. W. M. van der Geld, and J. G. M. Kuerten. Determination of the coefficients of langevinmodels for inhomogeneous turbulent flows by three-dimensional particle tracking velocimetry and direct numericalsimulation.
Physics of Fluids , 19, 2007.[39] A. Arn`eodo, R. Benzi, J. Berg, L. Biferale, E. Bodenschatz, A. Busse, E. Calzavarini, B. Castaing, M. Cencini,L. Chevillard, R. T. Fisher, R. Grauer, H. Homann, D. Lamb, A. S. Lanotte, E. L´ev`eque, B. L¨uthi, J. Mann,N. Mordant, W.-C. M¨uller, S. Ott, N. T. Ouellette, J.-F. Pinton, S. B. Pope, S. G. Roux, F. Toschi, H. Xu, andP. K. Yeung. Universal intermittent properties of particle trajectories in highly turbulent flows.
Phys. Rev. Lett. ,100:254504, Jun 2008.[40] F. Toschi and E. Bodenschatz. Lagrangian properties of particles in turbulence.
Annual Review of Fluid Me-chanics , 41(1):375 – 404, 2009.[41] Alex Liberzon, Beat L¨uthi, Markus Holzner, Søren Ott, Jacob Berg, and Jakob Mann. On the structure ofacceleration in turbulence.
Physica D: Nonlinear Phenomena , 241(3):208 – 215, 2012. Special Issue on SmallScale Turbulence.[42] R. Scatamacchia, L. Biferale, and F. Toschi. Extreme events in the dispersions of two neighboring particles underthe influence of fluid turbulence.
Phys. Rev. Lett. , 109:144501, Oct 2012.[43] A. Di Bernardino, P. Monti, G. Leuzzi, and G. Querzoli. Water-channel estimation of eulerian and lagrangiantime scales of the turbulence in idealized two-dimensional urban canopies.
Boundary-Layer Meteorology , 2017.[44] N. Stelzenmuller, J. I. Polanco, L. Vignal, I. Vinkovic, and N. Mordant. Lagrangian acceleration statistics in aturbulent channel flow.
Physical Review Fluids , 2:054602, May 2017.[45] J.I. Polanco, I. Vinkovic, N. Stelzenmuller, N. Mordant, and M. Bourgoin. Relative dispersion of particle pairsin turbulent channel flow.
International Journal of Heat and Fluid Flow , 71:231 – 245, 2018.[46] Ron Shnapp and Alex Liberzon. Generalization of turbulent pair dispersion to large initial separations.
Phys.Rev. Lett. , 120:244502, Jun 2018.[47] A. Celani, M. Cencini, M. Vergassola, E. Villermaux, and D. Vincenzi. Shear effects on passive scalar spectra.
Journal of Fluid Mechanics , 523:99108, 2005.[48] E. Pitton, C. Marchioli, V. Lavezzo, A. Soldati, and F. Toschi. Anisotropy in pair dispersion of inertial particlesin turbulent channel flow.
Physics of Fluids , 24(7):073305, 2012.[49] R. Shnapp, E. Shapira, D. Peri, Y. Bohbot-Raviv, E. Fattal, and A. Liberzon. Extended 3D-PTV for directmeasurements of Lagrangian statistics of canopy turbulence in a wind tunnel.
Scientific Reports , 9(7405), 2019.[50] Y. Bohbot-Raviv, R. Shnapp, A. Liberzon, V. Babin, M. Hotoveli, A. Shick, and E. Fattal. Turbulence statisticsof canopy-flows using novel lagrangian measurements within an environmental wind tunnel.
Physmod, LHEEA-DAUC, Ecole Centrale de Nantes , 2017.[51] C. W. Gardiner. Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences, 1997.[52] J. L. Doob. The brownian movement and stochastic equations.
Annals of Mathematics , 43(2):351–369, 1942.[53] M. Obukhov, A. Description of turbulence in terms of lagrangian variables. volume 6 of
Advances in Geophysics ,pages 113–116. Elsevier, 1959.[54] A. S. Monin and A. M. Yaglom.
Statistical Fluid Mechanics . Dover Publications inc., 1972.[55] M. S. Borgas and B. L. Sawford. Stochastic equations with multifractal random increments for modeling turbulentdispersion.
Physics of Fluids , 6(2):618–633, 1994.[56] J. D. Wilson, E. Yee, N. Ek, and R. d’Amours. Lagrangian simulation of wind transport in the urban environment.
Quarterly Journal of the Royal Meteorological Society , 135(643):1586–1602, 2009.[57] T. Duman, G. G. Katul, M. B. Siqueira, and M. Cassiani. A velocity–dissipation lagrangian stochastic model forturbulent dispersion in atmospheric boundary-layer and canopy flows.
Boundary-Layer Meteorology , 152(1):1–18,Jul 2014.[58] B. L. Sawford. Reynolds number effects in lagrangian stochastic models of turbulent dispersion.
Physics of FluidsA: Fluid Dynamics , 3:1577, 1991. [59] S. Du, B. L. Sawford, J. D. Wilson, and D. J. Wilson. Estimation of the kolmogorov constant (c0) for thelagrangian structure function, using a secondorder lagrangian model of grid turbulence. Physics of Fluids ,7(12):3083–3090, 1995.[60] Th. Dracos.
Three-Dimensional Velocity and Vorticity Measuring and Image Analysis Technique: Lecture Notesfrom the short course held in Zurich, Switzerland . Kluwer Academic Publisher, 1996.[61] OpenPTV consortium. Open source particle tracking velocimetry, 2014.[62] Y. Meller and A. Liberzon. Particle data management software for 3dparticle tracking velocimetry and relatedapplications – the flowtracks package.
Journal of Open Research Software , 2016.[63] R. Shnapp. Wind tunel canopy flow 3d-ptv. YouTube, 2018. .[64] A. N. Kolmogorov. The local structure of turbulence in incompressible viscous fluid for very large reynoldsnumbers.
Cr Acad. Sci. URSS , 30:301–305, 1941.[65] K. R. Sreenivasan. On the universality of the kolmogorov constant.
Physics of Fluids , 7:2778 – 2784, 1995.[66] M. Chamecki and N. L. Dias. The local isotropy hypothesis and the turbulent kinetic energy dissipation rate inthe atmospheric surface layer.
Quart. J. Roy. Meteor. Soc. , 130:2733–2752, 2004.[67] D. Poggi and G. G. Katul. Evaluation of the turbulent kinetic energy dissipation rate inside canopies by zero-and level-crossing density methods.
Boundary-Layer Meteorology , 136(2):219–233, Aug 2010.[68] S. B. Pope.
Turbulent Flows . Cambridge University Press, 2000.[69] H. Tennekes and J. L. Lumley.
A First Course in Turbulence . The MIT Press, 1972.[70] N. Mordant, E. L´evˆeque, and J.-F. Pinton. Experimental and numerical study of the lagrangian dynamics ofhigh reynolds turbulence.
New Journal of Physics , 6(116), 2004.[71] L. Biferale, E. Bodenschatz, M. Cencini, A. S. Lanotte, N. T. Ouellette, F. Toschi, and H. Xu. Lagrangianstructure functions in turbulence: A quantitative comparison between experiment and direct numerical simulation.
Physics of Fluids , 20(6):065103, 2008.[72] W.K. George, P.D. Beuther, and J.L. Lumley. Processing of random signals. In
Proceedings of the Dynamic FlowConference 1978 on Dynamic Measurements in Unsteady Flows . Springer, Dordrecht, 1978.[73] G. I. Taylor. Diffusion by continuous movements.
Proceedings of the London Mathematical Society , 1921.[74] C. Vanderwel and B. Ganapathisubramani. Turbulent boundary layers over multiscale rough patches.
Boundary-Layer Meteorology , 172(1):1–16, Jul 2019.[75] R. H. Shaw, Y. Brunet, J. J. Finnigan, and M. R. Raupach. A wind tunnel study of air flow in waving wheat:Two-point velocity statistics.
Boundary-Layer Meteorology , 76(4):349–376, Dec 1995.[76] P. K. Yeung and S. B. Pope. Lagrangian statistics from direct numerical simulations of isotropic turbulence.
Journal of Fluid Mechanics , 207:531–586, 1989.[77] N. Mordant, P. Metz, O. Michel, and J.F. Pinton. Measurement of lagrangian velocity in fully developed turbu-lence.
Physical Review Letters , 87, 2001.[78] B. L. Sawford and P. K. Yeung. Kolmogorov similarity scaling for one-particle lagrangian statistics.
Physics ofFluids , 23(9):091704, 2011.[79] H. Yu, K. Kanov, E. Perlman, J. Graham, E. Frederix, R. Burns, A. Szalay, G. Eyink, and C. Meneveau. Studyinglagrangian dynamics of turbulence using on-demand fluid particle tracking in a public turbulence database.
Journal of Turbulence , 2012.[80] R.-C. Lien and E. A. DAsaro. The kolmogorov constant for the lagrangian velocity spectrum and structurefunction.
Physics of Fluids , 14(12):4456–4459, 2002.[81] B. L. Sawford, P. K. Yeung, and J. F. Hackl. Reynolds number dependence of relative dispersion statistics inisotropic turbulence.
Physics of Fluids , 20(6):065111, 2008.[82] Y. Li, E. Perlman, M. Wan, Y. Yang, R. Burns, C. Meneveau, R. Burns, S. Chen, A. Szalay, and G. Eyink.A public turbulence database cluster and applications to study lagrangian evolution of velocity increments inturbulence.
Journal of Turbulence , 2008.[83] D. Poggi, A. Porporato, L. Ridolfi, J. D. Albertson, and G. G. Katul. The effect of vegetation density on canopysub-layer turbulence.
Boundary-Layer Meteorology , 111(3):565–587, Jun 2004.[84] H. Xia, N. Francois, H. Punzmann, and M. Shats. Lagrangian scale of particle dispersion in turbulence.
NatureCommunications , 4, 06 2013.
APPENDIX A - SUB-VOLUME FLOW PARAMETERS
The two following tables present the values of the sub-volume parameters that were estimated accordingto the description in Section II D. The change of Re λ with height is also presented in Fig. 8.9 sv ˜ u [m/s] (cid:15) [W/kg] η [mm] τ η [s] λ [mm] Re λ H/ηa a a a a b b b b b c c c c c d d d d d TABLE I: Turbulence parameters for each sub-volume for the Re ∞ = 16 × case. sv ˜ u [m/s] (cid:15) [W/kg] η [mm] τ η [s] λ [mm] Re λ H/ηa a a a a b b b b b c c c c c d d d d d TABLE II: Turbulence parameters for each sub-volume for the Re ∞ = 26 × case.FIG. 8: Change of the Taylor microscale Reynolds number with height for the various sub-volumes. fullsymbols are for Re ∞ = 1 . × , and open symbols are for Re ∞ = 2 . × .0 APPENDIX B - EMPIRICAL ESTIMATION OF AUTOCORRELATION FUNCTIONS
The autocorrelation of the random signals in this work were calculated as follows. Consider the set of i = 1 . . . N random series samples a i ( τ ) of the random variable a ( τ ). Generally speaking, the average of a ( τ ) and its standard deviation may change with τ , where τ = t − t , and t is the time at which the recordof a i began). The average of a ( τ ) is defined as µ ( τ ) = 1 N N (cid:88) i =1 a i ( τ ) (19)the fluctuations relative to the average are denoted a i ( τ ) (cid:48) = a i ( τ ) − µ ( τ ), and the standard deviation of a ( τ )is defined σ ( τ ) = (cid:34) N N (cid:88) i =1 a (cid:48) i ( τ ) (cid:35) / . (20)Note that these two definitions correspond to the sub-volume average introduced in section II C. Then theautocorrelation of a is calculated as follows ρ ( τ ) = N (cid:80) Ni =1 [ a (cid:48) i (0) a (cid:48) i ( τ )] σ (0) σ ( τ ) (21)The estimator Eq. (D.3) uses the average and the standard deviations that are allowed to change with τ . As discussed by Guala et al. [37], the Lagrangian trajectories with long tracking time possibly belong toa subset of “weak turbulence”. Consequently, an estimator of ρ that uses a single value µ and σ averagedover all values of τ is a biased estimator that may under predicts ρ ( τ ) at long times. Therefore using thedefinition σ ( τ ), that changes with τ in Eq.(D.3), prevents this underestimation at long time lags. This issuewas discussed in details by Guala et al. [37], where the biased estimator in their paper was denoted Eq. (2.1),and the unbiased estimator Eq. (D.3) here is equivalent to their Eq.(2.6).The autocorrelation functions were calculated in this article using many samples that were measuredduring the long experimental runs we have conducted, ∼ − −
15 minutes each. To demonstrate that ourestimations of the autocorrelation function and the decorrelation timescale are converged we show in themain panel of Fig. 9 the Lagrangian autocorrelation function with error bars that represent the results ofa bootstrapping calculation. Specifically, the dataset of trajectories in sub-volume b3 were divided to threegroups and the autocorrelation ρ xx was calculated separately. The error bars show the range of scatter ofthe results for the three groups and represents a small degree of uncertainty in the range relevant for ourstudy. Furthermore, the inset shows the convergence of T i that was calculated using subsamples of our datawith different sizes. The relative error of T i is seen to decrease rapidly with the subsample size. Therefore,Fig. 9 demonstrates that the autocorrelation and the decorrelation times were converged in our experiment,and suggests an uncertainty of up to a few percents.1 .
000 0 .
025 0 .
050 0 .
075 0 . τ [ms]0 . . . . . ρ x ( τ ) . . . . . . c o n v e r g e n ce o f T i FIG. 9: Convergence plots for the autocorrelation function (main panel) and for the decorrelation timescale(inset) for trajectories in sub-volume b3. The error bars in the main panel represent the range of scatterfor ρ xx calculated using 3 subsamples of the data. The inset shows the relative error in estimating T i usingsubsamples of different sizes. APPENDIX C - FINITE VOLUME EFFECT ON CORRELATION
Since the volume of observation is finite and due to the fact that occupation times of particles within thefinite volumes are dependent on their velocity, a natural bias occurs in the estimation of Lagrangian velocityautocorrelation functions in PTV experiments. To minimize this effect on the results shown in this work, theestimation of Lagrangian timescales was performed in this work only on short times, such that most of theparticles do not have sufficient time to leave the observation volume. A time scale for the occupation timeswithin a volume of dimension L is T vol = Lu (cid:48) (22)with u (cid:48) being the root mean squared value of particle velocities. In Fig. 10, the Lagrangian autocorrelationfunction for the x velocity component is presented against time normalized by T vol for the HIT DENS dataover two ensembles. The first is the full series of velocities over all trajectories. The second ensemble wasobtained by truncating the velocity series of each trajectory such that only values measured within a certainvolume of size L were taken to mimic the finite volume effect. The figure shows that for times in the range τ < T vol , the difference that exists between the two autocorrelation estimations is rather small - up to ∼ τ < . T vol , with L = 3mm and u (cid:48) determined from all samples at a given sub-volume. / T vol () full volumefinite volume FIG. 10: Lagrangian velocity autocorrelation over the full data set compared with trajectories that weretruncated to be within finite volume. Plotted against time normalized by the volume timescale, Eq. (22).2
APPENDIX D - ESTIMATION OF 3D EULERIAN STATISTICS
Using our Lagrangian dataset, we estimated Eulerian velocity statistics, such as the mean velocity andthe turbulent stresses by three-dimensional interpolations of the sub-volume averaged data. The scheme wasperformed as follows:1. we estimated a sub-volume averaged value for each desired statistics (e.g. (cid:104) u i (cid:105) = (cid:104) v i (cid:105) , and R ij = (cid:104) v (cid:48) i v (cid:48) j (cid:105) ).2. we obtained vertical profiles of each statistics above each sub-volume group (a, b, c, d, see Fig. 1(a))by using a linear interpolation with respect to z .3. we obtain 3D field estimates by using a two-dimensional inverse distance weighted interpolation of thefour vertical profiles.As an example, we present a two dimensional cut of the mean velocity field in Fig. 11. .
25 0 .
50 0 .
75 1 .
00 1 . x/H . . . . . z / H y = − . H . . . . . . | V | [ m / s ] FIG. 11: Two dimensional projection of the mean velocity over a plane parallel to the ( x, zx, z