Do coherent structures organize scalar mixing in a turbulent boundary layer?

Abstract A scalar emanating from a point source in a turbulent boundary layer does not mix homogeneously, but is organized in large regions with little variation of the concentration: uniform concentration zones. We measure scalar concentration using laser-induced fluorescence and, simultaneously, the three-dimensional velocity field using tomographic particle image velocimetry in a water tunnel boundary layer. We identify uniform concentration zones using both a simple histogram technique, and more advanced cluster analysis. From the complete information on the turbulent velocity field, we compute two candidate velocity structures that may form the boundaries between two uniform concentration zones. One of these structures is related to the rate of point separation along Lagrangian trajectories and the other one involves the magnitude of strong shear in snapshots of the velocity field. Therefore, the first method allows for the history of the flow field to be monitored, while the second method only looks at a snapshot. The separation of fluid parcels in time was measured in two ways: the exponential growth of the separation as time progresses (related to finite-time Lyapunov exponents and unstable manifolds in the theory of dynamical systems), and the exponential growth as time moves backward (stable manifolds). Of these two, a correlation with the edges of uniform concentration zones was found for the past Lyapunov field but not with the time-forward future field. The magnitude of the correlation is comparable to that of the regions of strong shear in the instantaneous velocity field.

The average p.d.f. of local scalar concentration in turbulence is a smooth near-Gaussian function, but our interest is in snapshots of the concentration field with size comparable to the boundary layer thickness. The associated p.d.f.s are marred by local maxima that signify large regions with near-homogeneous concentration. Identifying these regions starts with these p.d.f.s. We discuss two methods, of which cluster analysis is the preferred one.
In this paper, we also explore the relevance of two different types of flow structures for the organization of the concentration field. The first one is related to the rate at which two close points separate in the velocity field. It involves an integration over a time interval T using the velocity gradient field along a Lagrangian trajectory. A scalar field is constructed from the logarithm of the spreading rate: the FTLE Λ T (x). Spreading can also be tracked for backward time, then, the FTLE Λ −T (x) measures the rate of convergence of two fluid parcels towards x. Ridges of local maxima of these scalar fields are known as Lagrangian coherent structures (LCS) (Shadden, Lekien & Marsden 2005). In 2-D time-dependent chaotic flows, LCS are analogous to separatrices, in time-in dependent chaotic fields, where they form impenetrable barriers for scalar transport. In time-dependent two-dimensional flows, LCS form transient barriers of transport, while they are advected almost passively. These objects offer the tantalizing prospect of predicting the spread of contaminants from knowledge of large-scale coherent structures only.
This concept has been applied in 2-D turbulence (Haller & Yuan 2000;Haller 2015) and effectively 2-D geophysical flows; see Bettencourt, Lopez & Hernandez-Garcia (2012) and references therein. Twardos et al. (2008) analysed an experiment on a weakly turbulent 2-D flow, using the related concept of stretching fields (Haller 2015). Line-shaped maxima of the finite-time past Lyapunov field were found to act as barriers of scalar transport. The finite-time Lyapunov field was measured in a turbulent boundary layer in a plane perpendicular to the boundary by Chong, Wang & Zhang (2009). The relatively small main velocity U ∞ ≈ 9 × 10 −2 m s −1 and small Re θ = 481 allowed the tracing of structures up to 2 s, and allowed identification of hairpin-like structures in Λ −T (x). Wilson, Tutkun & Cal (2013) identify LCS in a plane parallel to the boundary from the 2-D velocity field measured in a turbulent boundary layer. By positioning the plane close to the boundary, the relevant integration time T is the local eddy turnover time, which could be met in the Eulerian frame. The observation time was further extended using Taylor's frozen turbulence hypothesis.
LCS may be a useful concept in fully developed 3-D turbulence. When they organize as sheets, LCS structures may hinder scalar mixing and may thus be associated with boundaries of UCZs. In this paper, we will explore this concept in an experiment of a turbulent boundary layer that provides the full 3-D velocity field. Using this information, we measure the finite-time Lyapunov field in a plane perpendicular to the boundary. In addition, we identify the edges of UCZs, and perform a conditional average of the finite-time Lyapunov field on these edges. The other flow structures studied in this paper are regions of strong shear. These zones organize momentum transport, and correlate with the boundaries of uniform momentum zones (Eisma et al. 2015). The question is whether these regions are also correlated with the edges of the UCZs.
LCS became popular in the context of 2-D (geophysical) flows. The question is how to fruitfully apply this interesting concept to scalar dispersion in fully developed 3-D turbulence.
We will first describe our experimental methods in § 2. Next, we explain in § 3 how we extracted UCZs, while the computation of finite-time Lyapunov fields and shear regions  are explained in § § 4.1 and 4.2, respectively. With these preliminaries, we then describe in § 5 the relation between coherent structures of the velocity field and the edges of UCZs.

Experiment
In our experiment we inject dye from a point into a turbulent boundary layer in a water channel. The set-up is illustrated in figure 1. The channel has a 5 m long experiment section and cross-section of 0.6 × 0.6 m 2 . Fluorescent dye (Rhodamine B) is injected from a point at the boundary, 0.75 m upstream from the measurement volume, with injection velocity equal to the local mean velocity, 0.5 U ∞ , where U ∞ is the free-stream velocity. This ensured minimal perturbation of the velocity field. The measurement volume, located 3.5 m from the boundary layer tripping point, is at one of the sidewalls. The 3-D velocity field in a box L x × L y × L z = 5.72 × 6.29 × 0.57 cm 3 was measured using tomographic PIV. In terms of the boundary layer height δ 99 = 3.81 cm the measurement volume is 1.5 × 1.65 × 0.15 δ 3 99 . Simultaneously, the 3-D concentration field was measured in five slices spanning the L z = 0.15 δ 99 spanwise side of the measurement volume. Both fields were registered using a scanning laser sheet at λ = 527 nm.
For the simultaneous measurement of the velocity and scalar concentration a scanning tomographic PIV system was combined with scanning LIF. Four high-speed CMOS cameras (Dimax, PCO) are used for tomographic PIV. Two cameras were equipped with f = 200 mm objectives, whereas the others were equipped with f = 105 mm objectives, all operating at f # = 11. This optics involved Scheimpflug adapters to allow for large viewing angles while keeping the particle images in focus. The top two PIV cameras were positioned at an angle of β = 30 • with respect to the normal of the light sheet, whereas both side-viewing PIV cameras were looking at an angle of β = 45 • . In order to minimize astigmatic aberrations, water-filled prisms were employed for all PIV cameras. LIF is performed using a single high-speed camera (Fastcam, Photron) equipped with a f = 105 mm objective at f # = 8.
Illumination was provided using a 50 mJ double pulsed Nd:YLF laser (Darwin Duo 80M, Quantronix). The laser beam was focused with a spherical lens (f = 1500 mm), and a 90 mm wide and 1.4 mm thick light sheet was subsequently formed by two cylindrical lenses with focal lengths of −25 and 90 mm. This light sheet spans the streamwise (x) and the wall-normal direction (y) of the turbulent boundary layer, while the scanning is performed in the spanwise direction (z). The measurement volume is scanned with five thin laser sheets of 1.4 mm thickness each, with an overlap of approximately 30 %. A scanning frequency f s of 640 Hz results in a recording frequency of the full measurement volume of 128 Hz.
The flow is seeded with 10 µm diameter tracer particles (Sphericell) which are nearly neutrally buoyant. The seeding density in the tomographic particle image velocimetry (TPIV) images is around 0.025 particles per pixel. The excitation of the used dye (Rhodamine B) matches the 527 nm operating wavelength of the laser. Separation of the PIV and LIF signals is obtained by employing lowpass filters (PIV cameras) and a highpass filter (LIF camera).
Calibration of the distorted LIF and PIV images was performed using a 3-D calibration plate (Type 11, LaVision). The mapping of the distorted images was done using third-order polynomials in the x-and y-directions, whereas linear interpolation is used in the z-direction. Self-calibration, as proposed by Wieneke (2005), was performed to enhance the accuracy of the particle reconstruction. After several refinements, the remaining disparities were typically of the order of 0.1-0.2 pixels.
The particle volumes were reconstructed from the PIV images using five iterations of the fast-MART algorithm (Elsinga et al. 2006). The particle volumes were 3-D cross-correlated using an iterative multi-grid scheme reaching a final interrogation box size of 32 3 voxels. This yields a spatial resolution of approximately 1.1 mm, corresponding to 34 wall units in each direction. From the relation for the Kolmogorov length η, η + = (κy + ) 1/4 , with κ = 0.41 (Stanislas, Perret & Foucaut 2008) we estimate η ≈ 0.2 mm at our largest distance from the wall. Our velocity resolution is ≈ 5 η; resolving η is probably not possible in TPIV. However, our interest is in large-scale structures. Velocity gradients were obtained by employing a second-order spatial regression filter ) with a filter size of 1.5 times the dimension of the interrogation box in each direction.
Care was taken to ensure that the concentration in the LIF images is sufficiently low for the fluid to be optically thin (Crimaldi 2008). In order to convert the local light intensity to a local dye concentration, a calibration procedure was performed. For this purpose, a small container with a known uniform concentration is positioned inside the measurement domain, after which images were recorded. This procedure was repeated for a few different known concentrations, resulting in a calibration curve for each pixel. Linear interpolation was employed to obtain the dye concentration from intensities at each pixel location. After calibration small negative concentrations were present in the free-stream region of the flow, which is a result of the finite accuracy of the LIF calibration. In order to remove these negative concentrations, the mean value of the concentration in the free-stream part of the flow is set to zero by subtracting it from the concentration field.
The turbulent characteristics of the boundary layer were measured with the TPIV set-up, with results shown in figure 2(a,b). Stereo PIV at a sampling rate of 40 Hz and collecting a total of 6 × 10 3 samples was used to check the results. Good agreement was found between the two methods. As figure 2 illustrates, the tomographic PIV results also agreed well with the data from DeGraaf & Eaton (2000).
At the measurement location the main stream velocity is U ∞ = 0.74 m s −1 , the boundary layer thickness is δ 99 = 38 mm. The resulting momentum and displacement thicknesses are θ = 4.0 mm and δ * = 5.2 mm, respectively, and the friction velocity is u τ = 31 mm s −1 . The momentum thickness Reynolds number is Re θ = 3050. The mean scalar concentrationc profile across the boundary layer is shown in figure 2(c). The maximum concentration is found around y/δ 99 = 0.1. The concentration approaches zero in the free-stream part of the flow, as expected. The maximum value of the mean concentration profile c m = 0.074 mg l −1 .
The 3-D velocity field in the L x × L y × L z = 1.5 × 1.65 × 0.15 δ 3 99 measurement volume was stitched together from the TPIV data in each of the 5 z-positions of the 1.4 mm thick laser sheet. The velocity field was corrected for the advection of the mean flow in the time interval (≈ 1.5 ms) between two subsequent laser sheet positions. The quality of the measured velocity field was assessed by computing its divergence. The true velocity field has zero divergence; for the experimental data we compute the normalized divergence computed over the entire velocity field which has zero mean and mean square variation ζ 2 1/2 = 0.21. The value of this quality measure is comparable to that in Casey, Sakakibara & Thoroddsen (2013), and references therein.
The z-resolution of the LIF images is determined by the width of the scanned laser sheet. In the present paper we focus on the concentration field in the central z = 0 plane. Therefore, this concentration field is an average over the 1.4 mm wide laser sheet. Much as for the scanned velocity field, also the concentration field has been corrected for the streamwise advection by the mean flow during a scan cycle.
Summarizing, the outcome of our experiment is a completely resolved 3-D velocity field together with a concentration field that is approximately resolved in slices of the 3-D measurement volume. This information is illustrated in figure 3. In the following sections we will quantify the relation between coherent structures of the concentration field and those of the velocity field. Figure 4 shows three snapshots of the concentration field, taken one large-eddy turnover time δ 99 /U ∞ apart. In these snapshots, UCZs can be recognized as large regions where the variation of the concentration is small, bordered by sharp gradients. The concentration  fields are from a 1.4 mm thick slice, centred in the spanwise extent of the 0.57 cm thick measurement volume. The idea of UCZs is intuitively appealing, but methods to find them must be well defined. Here, we follow the well-established methods of cluster analysis (Bezdek 1980), which were recently introduced in turbulence (Fan et al. 2019). We will compare our results with the method used by de Silva, Hutchins & Marusic (2016), Adrian, Meinhart & Tomkins (2000) and Eisma et al. (2015).

Uniform concentration zones
The starting point is the histogram P his (c) of concentrations. In the case of n zones clustered concentration values (which may still be spread over the observation plane), this histogram is characterized by n zones well-separated peaks. Each cluster corresponds to a peaked distribution. The concentrations in each of those still vary, but the variation is much smaller than that in the entire snapshot. In practice, clusters are not so well defined and the art becomes estimating the number n zones of clusters and to decompose the histogram into a sum of peaked distributions that are centred on n zones distinct concentration values. The first task is done using the statistical technique of kernel density estimation (Silverman 1986), while the tessellation of the concentration field in zones is done using the fuzzy cluster method (Bezdek 1980).
The histogram P his (c) is a discrete estimation of the underlying p.d.f. P(c). The number of clusters is the number of local maxima of P(c). A standard statistical procedure exists to overcome ambiguity associated with the width and placement of the discrete histogram bins (Silverman 1986). The trick is to endow each pixel value with a Gaussian distribution of concentrations with width h, and sum them. A question is the choice of the smoothing factor h. It can be proven that, for a Gaussian probability distribution function P(c), the optimum h is h ≈ σ P n −1/5 , where σ P is the standard deviation of P(c). Of course, in our case P(c) is not a Gaussian (we are counting its local maxima), but experience teaches us that this choice for h is adequate.
Cluster analysis to find UCZs can be compared with the ad hoc method for finding uniform momentum zones described by de Silva et al. (2016), Adrian et al. (2000) and Eisma et al. (2015). Let us now briefly describe this method. For each concentration snapshot c(x, t) a histogram P his (c i ) is constructed of n his = 32 equally spaced concentrations c i , i = 1, . . . , n his that span the dynamic range of c(x, t). The number n his is chosen small enough to counter noise, but large enough to distinguish peaks. Then, (i) the maximum c max , while its width σ is determined in a least squares procedure. (iii) This peak is deleted from the histogram: P his (c i ) → max(P his (c i ) − P G (c i ), 0), after which steps (i) through (iii) are repeated n zones times. (iv) After sorting the Gaussians with respect to their peak position, the intersections of neighbouring Gaussians are found. These intersections constitute the boundaries of the UCZs. Notice that this assignment of zones differs from that in the cluster analysis: concentration in the tail of a Gaussian is now associated with another cluster. No such ambiguity exists in the cluster analysis, which is also insensitive to the choice of discrete histogram bins. The two methods described will find different zone edges, and will find a different association with flow structures.
The average number of zones from the kernel density estimate is n zones = 3.7 ± 0.9, with 0.9 the root-mean-square variation. Therefore, we fixed the number of zones to four. An estimate of the number of zones from the discrete histogram is fraught with uncertainty. We will always rank the zones according to their area. Thus, the 'blue sky' in figure 5 has rank 1, and the neighbouring (green) zone has rank 3.
The outcome of both methods for a selected concentration snapshot is illustrated in figure 5; both methods find approximately the same zones, but there are differences. Below we will argue that the cluster method is superior, as it can discriminate more acutely between flow structures.
Apart from the technique of kernel density estimation, we believe that there does not exist an objective way to estimate the number of uniform zones. Moreover, as has been argued by Fan et al. (2019) in the context of uniform momentum zones, adding or subtracting 1 from the number of zones did not change their identification of interfacial layers.

Coherent structures of the velocity field
4.1. Finite-time Lyapunov exponents FTLEs gauge the exponentially fast spreading of nearby points. They were computed from the measured velocity field u(x, t) and strain field A(x, t), (A = ∇u, components A i,j = ∂u i /∂x j ), by integrating the evolution of small separations δ(t) along a Lagrangian trajectory, The time integration of (4.1) over a time interval [t 0 , t] defines the evolution matrix . This matrix was computed on a 2-D grid of initial points x 0 = x(t 0 ) in the central (z = 0) plane of the measurement volume. As time progresses, fluid parcels may leave this plane and explore other z values; these parcels are tracked in our 3-D velocity field. Therefore, the full 3-D measured velocity field was used in the computation of Λ ±T . Choosing the central plane minimizes the loss of fluid parcels exiting in the spanwise direction. The resulting Lyapunov field at z = 0 is compared with the concentration field of the central slice; the information in the other four slices was not used. The Lyapunov exponent is related to the separation of infinitesimally close points. When a computation is done using actual points advected by the velocity field u(x, t), their separation has to remain small, which necessitates their replacement when their separation grows too large. If not, the advected points which started near x 0 may move to completely uncorrelated regions of the velocity field and no longer reflect the Lyapunov exponent for the Lagrangian trajectory that started at x 0 . This is irrelevant if instead the gradient field A is used, but the statistical noise of A is larger than that of u. The accuracy of measuring gradients of the velocity field in our experiment has been addressed through the divergence error (2.1). The integration of (4.1) was done using a simple forward Euler scheme, and the fields at x(t) were computed from the measured 3-D PIV velocity fields using trilinear interpolation.
The largest eigenvalue λ 3 of the positive Cauchy-Green tensor with † the adjoint, defines the finite-time future Lyapunov exponent Λ T as with T = t − t 0 . Our measurements of the velocity field are Eulerian, which means that the time interval T is restricted by the advection of fluid parcels until they exit the field of view. Consequently, there is a distribution of Lagrangian times T over the field of view. Most of the time, trajectories leave the downstream side of the field of view, but they may also exit a z-boundary of the measurement volume. In figure 6 we show these times for one snapshot of Λ −T . (Throughout we use the shorthand Λ T for Λ t 0 +T t 0 and Λ −T for Λ t 0 −T t 0 .) Backward in time moves us upstream, so that the integration times are shortest for the left part of the frame. With increasing integration times structures of Λ ±T become narrower.
The Lyapunov field Λ T (x) measures the spreading of nearby fluid parcels in the future. However, the boundaries of UCZs were formed in the past. The past of fluid structures can be studied by starting Lagrangian trajectories at x 0 in the field of view, and integrating equation (4.1) backward in time. The corresponding Lyapunov field Λ t 0 −T t 0 (x) is then again defined in terms of the eigenvalue λ 3 .
Note that the smallest eigenvalue of C t 0 −T t 0 is comparable to the forward-time Lyapunov exponent at time t 0 − T. Local maxima of this field indicate regions in the flow where points were separated in the past, but have converged to x 0 at the instant of observation. On the other hand, structures are advected, so that future Lyapunov exponents already existed in the past and then may have organized the concentration distribution. The question then is whether it is the past or future Lyapunov field that is most relevant for boundaries of the UCZs. Another question concerns the time t at which the scalar field c(x, t ) should be compared with structures of C t+T t or C t−T t : should t be chosen t or a later or earlier time?
When computing statistics, we restrict field coordinates (x, y) to the upstream half of the observation window for past Lyapunov exponents, and the downstream half of the future Lyapunov exponents. This corresponds to integration times |T| 0.05 s. In finding the correlation of Λ ±T with the edges of UCZs we only include points where Λ ±T is convex in the direction of the eigenvector ξ 3 corresponding to the largest eigenvalue λ 3 of C t±T t , ξ 3 · H · ξ 3 < 0, with the Hessian H ij = ∂ 2 Λ ±T (x)/∂x i ∂x j , and ξ 3 the projection of ξ 3 on the xy-plane. More steps to refine the FTLE field to ridges were mentioned by Farazmand & Haller (2012). We finally impose a threshold Λ ±T > 4 s −1 . This threshold corresponds to the noise level of Λ −T in the 'blue sky' outside the boundary layer (zone 1 in figure 5).

Shear vorticity ω sh
As a second velocity field candidate for the edges of UCZs, we consider regions of strong shear. In order to isolate those, the 2-D (u, v) projection of the velocity field was decomposed into three parts (Kolár 2007). This decomposition was introduced originally to provide a more robust description of vortices in turbulent flows, but is now used to detect shear layers (Maciel, Robitaille & Rahgozar 2012).
With this aim, the velocity gradient ∇u is separated into components representing rigid-body rotation, elongation and the desired true shear. The shear vorticity ω sh is then defined as (4.5) where u and x refer to a coordinate system formed by the principal axes of ∇u, rotated over 45 • . Local maxima of this field are the shear layers found by Meinhart & Adrian (1995), which separate regions with nearly uniform velocity (Eisma et al. 2015). A generalization of (4.5) does not exist, but an alternative way to find shear layers, but now using the full 3-D velocity field, has been proposed by Horiuti & Takagi (2005). We also would like to point to Haller & Beron-Vera (2012), who propose quantifying shear using the eigenvectors of the Cauchy-Green tensor (4.2). Summarizing, our measurement of the space-time velocity field in a 3-D volume allows us to compare three candidate velocity structures. The past Lyapunov field Λ −T (x), which relates to the past organization of the concentration field and attracting structures, the 2-D shear vorticity ω sh (x), which refers to the instantaneous velocity field and the future Lyapunov field Λ T (x) which corresponds to repelling structures. Figure 7 shows the concentration, the finite-time Lyapunov exponent, and the shear vorticity in a snapshot of the xy plane centred in the measurement volume. We register time series at a frequency of 128 Hz (≈ 6 U ∞ /δ 99 ), lasting 6.5 s, corresponding to 1.2 × 10 2 large-eddy turnover times δ 99 /U ∞ . In collecting statistics, a narrow region near the wall, where gradients are large, was excluded. This region is indicated in figure 7; its width is 150 y + (4.73 mm), which includes the peak in u in figure 2. Our experiments are in the Eulerian reference frame. Due to limitations on the integration time, the field Λ −T in the left side of figure 7(b) is more diffuse than in the right side of the figure.

Results
We will now discuss the relation between the Λ −T and ω sh fields, and finally present our main result, that is the relation between maxima of the finite-time Lyapunov field and edges of UCZs.

Statistics of Λ ±T and ω sh
Loosely speaking, Λ ±T is the time average of the largest stretching induced by the strain A over a Lagrangian trajectory, while ω sh is the instantaneous value of the strongest shear that it induces. The p.d.f.s of Λ −T and Λ T in figure 8(a) are remarkably similar, but that of ω sh has a long tail. Local maxima of all three fields are candidate structures organizing the distribution of passive scalar. Still, their nature is very different. Not surprising, the fields Λ −T and ω sh are not strongly correlated, as is illustrated in figure 8: for each ω sh , Λ −T has a broad distribution, with the average Λ −T increasing very slightly with increasing ω sh . The regions of large ω sh in figure 7 are diffuse blobs, so are those of Λ −T . Sharp line-like structures of Λ ±T are only expected for long enough integration times T, which should be of the order of the large-scale turnover time. On average, our maximum integration time T ≈ 0.1 s, while the typical eddy turnover time δ 99 /u rms ≈ 1 s. Clearly, sharply defined structure of Λ ±T needs a true Lagrangian measurement. Nevertheless, it is possible to measure the statistics of blobs with large Λ ±T and ω sh . By computing the second-order moment matrix of these regions and finding the eigenvector with the largest eigenvalue, the orientation of these structures can be found. As figure 9 illustrates, patches of large Λ ±T ≥ 20 s −1 are inclined at an angle of ϕ ≈ 20 • , which agrees with Chong et al. (2009) who identified these regions with hairpin vortices. Patches of large Λ T do not have a preferential orientation whereas the shear field ω sh appears to be oriented parallel to the boundary. As will be argued below, edges of concentration zones correspond to Λ −T , and much less to Λ T and the association of Λ −T (and not Λ T ) with genuine structures is tempting.

Conditional averages
In finding a connection between the edges of uniform scalar concentration zones and regions of large Λ ±T we follow the procedure of Bisset, Hunt & Rogers (2002), but with an important caveat. The procedure is to compute the conditionally averaged function Λ ±T (x, y − y 0 ; t) x,t , with y 0 the vertical (wall-normal) coordinate of an edge of a uniform concentration zone. The caveat is that we compute conditional averages relative to those where y 0 is chosen randomly. Specifically, we show the ratio of the average function Λ ±T (x, y − y 0 ; t) x,t , with y 0 the true zone edge over the average function with y 0 taken randomly. When this ratio is one, no distinction can be made between the true and random y 0 , and zone edges are not correlated with Λ ±T . The random choice of y 0 is such that its p.d.f. is the same as that of the true y 0 . We will denote the relative averages as Λ −T,T and ω sh . Further details are in the Appendix. Conditional averages are shown in figure 10. We have computed relative conditionally averaged vertical profiles of Λ ±T and ω sh on the boundaries of the second largest UCZ shown in figure 10(d), which corresponds the green zone in figures 5(b,c) and 7(a), but the conclusions are comparable to those involving all boundaries. In addition, we have removed small isolated patches (relative area < 5 × 10 −3 of the full frame) of this zone and zones of other rank. Finally, we recall that we only consider Λ −T for the right half of the frame, and vice versa for Λ T .
The central result of this paper, shown in figure 10, involves the conditional statistics of Λ ±T and ω sh on zone edges. Those edges are determined in two ways: using the cluster method (figure 10b,c,e, f,g) and using the simple histogram method (figure 10d). In the case of the cluster method, the correlation of the past Lyapunov field Λ −T with concentration edges is much larger than that of the future Λ T (figure 10b,c). This agrees with the result of Twardos et al. (2008) in weakly turbulent 2-D flow who found that contours of an advected scalar field follow the past stretching field, and that of Kushnir et al. (2006) in numerical turbulence (Re λ = 24) who conclude that gradient sheets of the dissipation field are formed by the contracting eigenvalue of M T . A similar correlation with edges of the concentration zones is seen for the relative shear vorticity ω sh in figure 10(e).
The increase at a zone edge of the past FTLE field Λ −T is a mere 0.05 (5 %) above the random background Λ −T = 1, and a similar amount above the future Λ T . The significance of this enhancement should be compared with the concentration of ω sh on edges of uniform momentum zones (regions of uniform u(x, t)), ω sh ≈ 1.3 (30 %), shown in figure 10(g). It was derived from the present data using the same cluster analysis. The linkage of these two structures -both derived from the velocity field -has been discussed in Eisma et al. (2015), Christensen & Adrian (2001) and Elsinga & Marusic (2010).
The significance of our result is further illustrated in figure 10(b,c) where the correlation is seen to depend sensitively on the time delay t between the FTLE field and the concentration zone edges. Specifically, we correlate edges of the concentration field c(x, t + t) with the FTLE fields Λ t+T t and Λ t−T t for t increasing from −0.023 s to 0 s. Although the field Λ t−T t explores the past, the conditional average peaks at the present. The cluster method needs the number of zones n zones as input. Throughout, we have fixed it to the averaged kernel density estimate, n zones = 4, both for the cluster and the histogram method. Figure 10 count the blue sky as a zone, and take conditional averages on the edges of the zone with rank 2. Therefore, fixing n zones = 3 almost always amounts to averaging with respect to the turbulent-non-turbulent interface. For the edges found with the histogram method (figure 10d) there is almost no distinction between past and future Lyapunov exponents; perhaps illustrating the superiority of the cluster method.
Summarizing, we find a small but significant correlation between the edges of UCZs, and both the past-time Lyapunov field and the shear vorticity field.

Conclusion
We have quantified the large-scale organization of scalar concentration in a turbulent boundary layer. In analogy with uniform momentum zones in turbulent boundary layers (Meinhart & Adrian 1995;Hutchins et al. 2012;Eisma et al. 2015) we have used the term 'uniform concentration zones', but we emphasize the strong connection with ramp-cliff structures (Warhaft 2000;Iyer et al. 2018;Buaria et al. 2021). Two different procedures were described for the identification of UCZs, with superior results being obtained by the cluster method. Still, it can miss zones, or can find zones which are not obvious to the eye. After all, the concentration inside a zone is not strictly uniform. Ramp-cliff structures share the sudden change of concentration (cliffs) with the edges of UCZs. However, cluster analysis cannot identify ramp-cliff structures.
Next, we considered structures of the velocity field as candidate zone edges. From the measured 3-D velocity field and its gradients we compute the intersection of the 3-D finite-time Lyapunov field Λ ±T with the central plane of the measurement volume. Due to the advection by the mean velocity, the longest time interval over which trajectories could be integrated is T ≈ 2 δ 99 /U ∞ . This may be too short to see a clear development of a Lagrangian coherent structure. Further refinement of the Λ −T field to its local maximum ridges may decrease the random background of the measured correlations.
Concentration edges found with the cluster method only correlate with the past finite-time Lyapunov field Λ −T . They select converging rather than diverging flow structures. No such distinction between past and future can be made for concentration edges found with the histogram method. Perhaps this points to the superiority of the cluster method.
The concentration of both Λ −T and ω sh fields on the edges of UCZs is quite small compared with the random background. Using our techniques, it is a mere 5 %. This value can be compared with the concentration of shear vorticity on edges of uniform momentum zones, which is well established in the literature; it is 30 % using the same techniques.
Our conclusion is that the question as to which large-scale velocity structures organize the scalar concentration is still open.  Figure 11. Conditional averages of finite-time Lyapunov exponents Λ ±T (x, y − y 0 ; t) x,t , and shear vorticity ω sh (x, y − y 0 ; t) x,t on both randomly picked y 0 and true vertical (wall-normal) coordinates y 0 of UCZs. The figure illustrates the necessity of relative conditional averages. (a) Compares conditional averages with respect on the true edges y 0 (black lines) with those on the randomly picked y 0 (grey lines). For ω sh , the displacement of this structure from y − y 0 = 0 is roughly the correlation length. (b) Profiles of Λ −T (x, y; t) x,t and Λ T (x, y; t) x,t (full lines) and ω sh (x, y; t) x,t (dashed line). (c) Profile P( y 0 ) of zone edges; it is obscured by its randomized version.

Appendix. Accidental conditional averages
We demonstrate that conditional averages may show a non-trivial dependence on y − y 0 , even for randomly picked y 0 . The reason is that the ordinary averages Λ ±T (x, y; t) x,t and ω sh (x, y; t) x,t strongly depend on the height y in a turbulent boundary layer. Conditional averages on the genuine edges y 0 should, therefore, be taken relative to ones using randomly picked y 0 . Because also the distribution of edges y 0 is inhomogeneous, the randomly picked coordinates y 0 should have the same p.d.f. as the genuine ones. This can trivially be achieved by using a weighted pseudo-random number generator. Although further refinements on this 'null hypothesis' are possible, we believe that this approach suffices to discriminate real from accidental correlations in our experiment. Figure 11(a) compares conditional averages with respect to real and random vertical coordinates of edges y 0 . Clearly, the random conditional averages show structure (steps), which is due both to the correlation properties of Λ ±T and ω sh , and to the non-homogeneous distribution of the averages Λ ±T (x, y; t) x,t and ω sh (x, y; t) x,t , which are shown in figure 11(b). Conditional averages with respect to the true edge locations y 0 involve the p.d.f. of y 0 , as do the averages with respect to the randomly picked y 0 . Both p.d.f.s are shown in figure 11(c). The p.d.f. of the edge locations is proportional to the mean concentration gradient |∂c(x, y; t)/∂y| x,t . To the best of our knowledge, the statistical bias of conditional averages seems to have been ignored in the literature so far (Bisset et al. 2002;de Silva et al. 2017).
The emergence of non-trivial conditional averages can be understood using a simple analytical argument. We demonstrate that a conditional average of random functions on random points shows structure, even though the points and functions are completely uncorrelated. This is the case when the average function itself has a non-trivial dependence on its argument.
Imagine we sample random functions f i ( y) at random points y 0 . In our case these random functions are Λ ±T (x, y) and ω sh (x, y) at different x and snapshots t, while they are conditioned on the vertical coordinate y 0 of edges.
We first do the average over i at fixed y 0 . Then, for a (small) interval around y 0 , we Taylor expand with respect to y = y − y 0 . The average function is f i ( y 0 ) + ( d/ dy 0 )f i ( y 0 ) y + 1 2 ( d 2 / dy 2 0 )f i ( y 0 ) y 2 + . . ., where denotes the average over realizations of the function f i ( y). Since differentiation is a linear operation, we may switch averaging and differentiation, so that the average over realizations is f i ( y 0 ) + ( d/ dy 0 ) f i ( y 0 ) y + 1 2 ( d 2 / dy 2 0 ) f i ( y 0 ) y 2 + . . .. We next average over randomly picked points y 0 . Now, the polynomial representation of the conditionally average involves ( y − y 0 ) n ( d n / dy n 0 ) f i ( y 0 ) y 0 . If the average function f i ( y 0 ) is not flat and has a few non-zero derivatives, the conditional average has structure.