2.1. Numerical models

LES solves the low-pass filtered Navier–Stokes equations for turbulent length scales larger than the characteristic grid dimension ${\it\Delta}$ . In the absence of non-inertial and buoyancy forces, the filtered Navier–Stokes equations can be expressed as (written in rotational form)

where a tilde ( ${\sim}$ ) denotes a filtered quantity, $\tilde{u} _{i}$ (using index notation) denotes the filtered velocity vector in each of the Cartesian directions $x_{j}$ ( $1=u=\text{streamwise}$ , $2=v=\text{spanwise}$ , $3=w=\text{vertical}$ ), $t$ is time, $\tilde{p}^{\ast }=\tilde{p}+(\tilde{u} _{j}\tilde{u} _{j})/2$ is the modified pressure, ${\it\nu}$ is the dynamic viscosity and $F_{i}$ is a generic forcing term. $D_{i}$ is the force resulting from the drag of (unresolved) canopy elements, which are aggregated over a cell volume. This is calculated according to a standard drag law (Shaw & Schumann Reference Shaw and Schumann1992)

where $c_{d}$ is the drag coefficient, $a$ is the one-sided leaf area density (LAD) and $\tilde{V}=(\tilde{u} _{i}\tilde{u} _{i})^{1/2}$ . The subgrid-scale (SGS) stress tensor, ${\it\tau}_{ij}$ , represents the effects of unresolved turbulent motions and must be modelled. In this study, this is done using the well-known Smagorinsky model, which calculates the SGS stress tensor as

where $\tilde{\unicode[STIX]{x1D61A}}_{ij}=((\partial \tilde{u} _{i}/\partial x_{j})+(\partial \tilde{u} _{j}/\partial x_{i}))/2$ is the resolved strain rate tensor and $|\tilde{\unicode[STIX]{x1D64E}}|=(2\tilde{\unicode[STIX]{x1D61A}}_{ij}\tilde{\unicode[STIX]{x1D61A}}_{ij})^{1/2}$ . $C_{s}$ is the Smagorinsky coefficient, which is scale dependent and calculated dynamically along fluid particle trajectories following Stoll & Porté-Agel (Reference Stoll and Porté-Agel2006a ). Test filtering for the scale-dependent scheme is performed at scales of $2{\it\Delta}$ and $4{\it\Delta}$ .

Passive scalar transport is calculated by solving the filtered scalar-transport equation, which is given by

where $\tilde{{\it\theta}}$ denotes a generic (filtered) scalar quantity, ${\it\alpha}_{{\it\theta}}$ is its molecular diffusivity, $S_{{\it\theta}}$ is a generic scalar source (flux divergence) and $q_{i}$ is the SGS scalar flux which is calculated as

where $Sc_{sgs}$ is the SGS Schmidt number, which is scale dependent and calculated dynamically following Stoll & Porté-Agel (Reference Stoll and Porté-Agel2006a ).

A numerical solution for (2.2) and (2.5) is obtained by discretizing the computational domain of size $L_{x}\times L_{y}\times L_{z}$ into uniform hexahedral elements of size ${\it\Delta}_{x}\times {\it\Delta}_{y}\times {\it\Delta}_{z}$ , which gives a characteristic grid scale of ${\it\Delta}=({\it\Delta}_{x}{\it\Delta}_{y}{\it\Delta}_{z})^{1/3}$ . Horizontal and vertical gradients are calculated using Fourier-based pseudospectral and second-order finite-difference schemes, respectively. Integration in time is performed using a second-order Adams–Bashforth scheme. Aliasing errors that result from nonlinear terms are removed using the $3/2$ rule (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang1988).

The computational domain (and therefore canopy) is periodic in the horizontal directions. The upper and lower boundaries are zero-flux rigid surfaces, where the upper boundary has zero stress and the lower boundary is a no-slip wall. The viscous stress at the lower wall is unresolved, and is therefore modelled by specifying the shear stress at the lowest computational level ( $z=0.5{\it\Delta}_{z}$ ) using Monin–Obukhov similarity theory (Stoll & Porté-Agel Reference Stoll and Porté-Agel2006b ). The flow is driven by a spatially-uniform streamwise pressure gradient. The passive scalar $\tilde{{\it\theta}}$ also follows a zero-flux condition at the wall. Different domain and discretization configurations were used for the start-up and quasi-steady simulations, which are described independently in §§ 3 and 5.

The LAD was horizontally homogeneous and was also chosen to be constant in the vertical to avoid possible complexities associated with profile shape such as the formation of a secondary wind-speed maximum in the canopy understorey. The LAD was reduced at the top two numerical grid points in the canopy to transition smoothly to zero at the canopy top, which was to minimize errors associated with the dispersive numerical scheme used to compute vertical gradients. For simplicity, the drag coefficient $c_{d}$ was held constant in space, although it may be more physically realistic to allow it to vary spatially (Finnigan & Shaw Reference Finnigan and Shaw2008) or as a function of wind speed or direction (Pan, Chamecki & Isard Reference Pan, Chamecki and Isard2014).

Simulation inputs for the start up and fully-developed cases were generally the same, with the primary difference being domain size. Specific input values are given separately in the following sections for each simulation case. Inputs were also chosen to be similar to that of Finnigan et al. (Reference Finnigan, Shaw and Patton2009) to facilitate direct comparison. The main differences between our simulations and that of Finnigan et al. (Reference Finnigan, Shaw and Patton2009) is the chosen SGS model, and that we use a larger domain height but a smaller horizontal domain extent. However, Finnigan et al. (Reference Finnigan, Shaw and Patton2009) noted little impact on results with increased horizontal domain size.

In what follows, an overbar denotes a time average over the longest available quasi-steady simulation period. A departure from this average is indicated by a prime, i.e. ${\it\phi}^{\prime }={\it\phi}-\overline{{\it\phi}}$ . Angle brackets $\langle \cdot \rangle$ denote an average at a given vertical level over the horizontal extent of the domain.

3.1. Roll up and vertical asymmetry

Previous authors have demonstrated that flows with an inflected velocity profile (e.g. mixing layers) are inviscidly unstable to disturbances of arbitrary wavenumber at all stages of development, with the most unstable being relatively low wavenumbers (Michalke Reference Michalke1965; Huerre & Monkewitz Reference Huerre and Monkewitz1985). As discussed earlier, the mean velocity profile of a canopy layer contains an inflection point at the canopy top which is a result of momentum absorption by canopy elements.

In the start-up canopy simulations (and analogous to mixing-layer flows) roll up occurs rapidly as a result of the growth of low-wavenumber disturbances. The canopy-top velocity inflection causes the formation of characteristic spanwise roller vortices, as illustrated by figure 2(a).

Symmetry is an important feature of mixing layers. The mixing layer grows symmetrically about the mixing plane, and velocity moments are generally symmetric (or antisymmetric) about the mixing plane (Bell & Mehta Reference Bell and Mehta1990; Rogers & Moser Reference Rogers and Moser1994). Previous work has demonstrated that the addition of a wall boundary to a mixing layer acts to break this symmetry (Lindzen & Rosenthal Reference Lindzen and Rosenthal1976; Holt Reference Holt1998). The addition of a boundary acts to displace the roller cores away from the mixing plane toward the unbounded side, although the roller still remains roughly symmetric (Holt Reference Holt1998). Figure 2(a) shows instantaneous contours of spanwise-averaged vorticity $\widetilde{{\it\omega}}_{2}=((\partial \tilde{u} /\partial z)-(\partial \tilde{w}/\partial x))$ for the start-up flow. In a free mixing layer, rollers are rotationally symmetric about the mixing plane. In contrast, rollers in the canopy layer are not symmetric about the canopy top (mixing plane). Figure 2(a) shows that the roller core is located at roughly $z=1.3h$ – $1.4h$ instead of $z=h$ as one might expect in symmetric flow. This symmetry breaking is to some degree a result of blocking by the lower wall as in Holt (Reference Holt1998). In the presence of a sufficiently dense canopy, the wall stress is essentially zero and the lower wall is not felt by upper canopy structures. In this case a displaced ‘effective’ wall is created at some height $0<z<h$ .

The canopy has an additional feature not present in previous studies of wall-bounded mixing layers, which is distributed momentum absorption by canopy elements. The canopy acts as a momentum sink and damps turbulent motions. Therefore, the portion of roller structures extending into the canopy experiences a non-conservative damping force that is not present above. In both cases, both mechanisms contribute to an asymmetrical eddy structure about the canopy top.

3.2. Scalar microfronts

As discussed previously, upward, ejection-generating coherent motions near the canopy top are inhibited by the low-flux canopy. However, downward, sweep-generating coherent motions are much stronger as they have a source of free air from above. Consequently, sweeps dominate ejections in the canopy. To satisfy mass or scalar continuity within the canopy, upward ejection events must be more frequent than the more powerful downward sweep events, which has been previously observed in experiments (e.g. Coceal et al. Reference Coceal, Dobre, Thomas and Belcher2007) and is demonstrated by our fully-developed LES (figure 7). In terms of the roller vortices, this means that a larger fraction of the canopy top surface area must be occupied by ejection events than sweep events.

The periodic cycle of sweeps and ejections gives rise to the commonly reported ‘scalar microfront’ structures, or what may be considered unusually large streamwise scalar gradients. Figure 2(e) shows instantaneous scalar perturbations along a streamwise transect at a height of $z/h=1$ for the simulated flow. The transect is a spanwise average taken after roughly 10 eddy turnover times from the beginning of the simulation. At this point in the flow evolution, three-dimensionality is very weak and the two-dimensional roller structure dominates the flow (see figure 4 for reference). Contours of the associated spanwise vorticity field are given in figure 2(a) which show the approximate location of the rollers. In the region extending from just upwind of the roller core to approximately midway between the rollers, transport is generally upward. Similarly, downward transport occurs in the region immediately downwind of the roller core.

A microfront occurs at the interface where transport abruptly transitions from downward to upward directions. It is important to note that the microfronts are not necessarily an indicator of a single locally-large vertical momentum flux. Rather they are a signature of opposing vertical fluxes in the upwind and downwind directions (figure 2 d). We confirm that the microfront coincides with local maxima in pressure (figure 2 b), as has been previously reported in other work examining fully-developed flow (Fitzmaurice et al. Reference Fitzmaurice, Shaw, Paw and Patton2004; Finnigan et al. Reference Finnigan, Shaw and Patton2009). Despite the fact that figure 2(e) indicates that the largest streamwise scalar gradients are between the rollers, figure 2(c) suggests that the largest (negative) turbulent stresses are generated in the sweep region upwind of the microfront and at the roller core.

What is additionally striking is the asymmetrical ‘ramp’-like patterns that form. Although scalar traces shown in figure 2(e) are associated with the two-dimensional, temporally-developing flow, they exhibit the classical features of experimental data collected in real canopy flows (e.g. Gao et al. Reference Gao, Shaw and Paw1989; Shaw et al. Reference Shaw, Paw and Gao1989). Figure 2(f) shows a sample scalar time series from the study of Gao et al. (Reference Gao, Shaw and Paw1989), which has a qualitatively similar trend as in the simulated flow. Both signals show a periodic trend of sharp streamwise gradients bracketed on either side by slowly decreasing scalar perturbations and a relatively quiescent zone. This offers compelling evidence that the scalar ramp patterns observed in real canopy flows could be caused by an underlying quasi-two-dimensional roller structure that persists to the fully-developed state.

In contrast to this result, previous studies of canopy structures have associated the microfronts with a more three-dimensional structure, such as the dual hairpin structure of Finnigan et al. (Reference Finnigan, Shaw and Patton2009). Such a three-dimensional structure appears to be supported by two-point correlation measurements performed by Shaw et al. (Reference Shaw, Brunet, Finnigan and Raupach1995), which found that the spanwise correlation scale is of order $h$ rather than the spanwise flow extent. However, it has been well documented that even slight irregularities in a predominantly spanwise structure can result in very low spanwise correlation scales. Browand & Troutt (Reference Browand and Troutt1980) and Browand & Troutt (Reference Browand and Troutt1985) demonstrated this fact for a predominantly two-dimensional mixing layer. Hutchins & Marusic (Reference Hutchins and Marusic2007) found that meandering of streamwise structures severely affected length scales inferred from two-point correlation measurements. We suspect that the near-canopy flow is dominated by a predominantly two-dimensional structure that contains irregularities. We certainly do not contend that the flow does not contain three-dimensionality. The presence of an approximate $k^{-5/3}$ kinetic energy spectrum observed in both fully-developed mixing and canopy layers (e.g. Huang & Ho Reference Huang and Ho1990; Finnigan Reference Finnigan2000) implies that the flow is fully three-dimensional. However, our claim is that underlying two-dimensional motions dominate near-canopy flow and are largely responsible for the scalar ramp structures. We now seek to more closely examine the role of three-dimensionality in the structure of the flow.

3.3. Development of three-dimensionality

As the start-up flow progresses, vortex breakdown occurs and smaller turbulent scales appear, as illustrated by the temporal progression shown in figure 4. At this point, the rollers are no longer well-defined vortex tubes as one may envision in the traditional sense. Rather, they are underlying vortical motions that transport smaller vortices in a circular motion. This process is similar to that observed in traditional mixing-layer studies (e.g. Tanahashi, Iwase & Miyauchi Reference Tanahashi, Iwase and Miyauchi2001; Wang, Tanahashi & Miyauchi Reference Wang, Tanahashi and Miyauchi2007). It is perhaps because of this fact that in flows with very high Reynolds numbers, Kelvin–Helmholtz rollers are usually referred to as ‘billows’ rather than vortices. Localized vortex sheets are generated by neighbouring structures, which can then undergo higher-order roll up, sustaining a self-similar structure (Moser & Rogers Reference Moser and Rogers1991).

Figure 4.

Progression of the start-up simulations as shown by isosurfaces of ${\it\lambda}_{2}$ . Isosurfaces are coloured by vertical height for visual clarity. Mean flow is from left to right. Panels show different stages of flow progression: (a) $t\,(U_{0}(h)/h)=1.23$ , ${\it\lambda}_{2}=-0.1$ ; (b) $t\,(U_{0}(h)/h)=16.8$ , ${\it\lambda}_{2}=-40$ ; (c) $t\,(U_{0}(h)/h)=18.5$ , ${\it\lambda}_{2}=-200$ . Frame (d) is the same as (c), except that isosurfaces are coloured by $\tilde{u} ^{\prime }/U(h)$ . An animation of the start-up simulation can be found in the online supplemental material (see supplementary movie 1 available at http://dx.doi.org/10.1017/jfm.2015.749).

As is required by the vorticity equation, vortex lines must continually stretch. Given sufficient time and a large enough Reynolds number, they are likely to form inclined arched or hairpin-like structures (Robinson Reference Robinson1991; Gerz et al. Reference Gerz, Howell and Mahrt1994). The precise inclination angle is flow-dependent and is the point when stretching, mean rotation and induced flow are in equilibrium (Rogers & Moin Reference Rogers and Moin1987). Gerz et al. (Reference Gerz, Howell and Mahrt1994) and Rogers & Moin (Reference Rogers and Moin1987) found in homogeneous turbulence and free shear flow, respectively, that hairpins tended to form in head-up and head-down vortex pairs. Accordingly, they also reported a relatively equal probability of detecting head-up and head-down hairpins. It is not surprising that they found equal probabilities of head-up and head-down hairpins given the symmetry of homogeneous and free shear flows.

As discussed previously, the wall breaks vertical symmetry and thus it would not be surprising to find preferential hairpin formation. In fact many previous workers have reported a more asymmetric hairpin vortex structure in the ISL, which leads to the dominance of head-up hairpin formation (cf. Adrian Reference Adrian, Venditti, Best, Church and Hardy2013).

Figure 4(c,d) shows that after the primary roller structures break down into smaller scales, the underlying two-dimensional structures remain. However, a highly three-dimensional structure also emerges as smaller vortices are stretched and strained, with the two evolving somewhat independently. This feature has been frequently reported in the late development of mixing layers (Browand & Troutt Reference Browand and Troutt1980, Reference Browand and Troutt1985; Lasheras & Choi Reference Lasheras and Choi1988). Superposed on the two-dimensional structure is an array of highly three-dimensional smaller-scale structures. Stretching combined with straining by the primary structures creates a convoluted entwinement of vortices. As is ubiquitously reported in mixing layers, initially spanwise structures stretch and are strained by the primary rollers to create arch- or hairpin-like structures oriented in the streamwise direction (Lasheras & Choi Reference Lasheras and Choi1988; Rogers & Moser Reference Rogers and Moser1992). Arched structures with their heads oriented downward and behind their legs are strained downward and against the mean streamwise flow. Their vertical progression is impeded by a combination of canopy drag and the wall (or displaced wall), but they still may form hairpins due primarily to horizontal straining. Head-up structures on the other hand are strained upward, and are not immediately bounded from above. Thus, we find that they are able to continue stretching upward into a classical boundary-layer hairpin, in many cases extending several canopy heights tall in the vertical.

We notice that the tallest head-up hairpins tend to form along common streamwise transects. For this particular flow, they are at $y=0=L_{y}$ and $y=L_{y}/2$ . If we colour the isosurfaces by streamwise velocity perturbation (figure 4 d), we find a distinctly ‘streaked’ pattern in the streamwise direction. The streaky patterns can also be confirmed by examining short time averages of the flow (not shown). This streamwise streaky structure that persists over very long streamwise distances has been well documented in mixing layers (Breidenthal Reference Breidenthal1978; Bernal et al. Reference Bernal, Breidenthal, Brown, Konrad and Roshko1979; Rogers & Moser Reference Rogers and Moser1992) and smooth-walled boundary layers (Kim & Adrian Reference Kim and Adrian1999; Jiménez, Álamo & Flores Reference Jiménez, Álamo and Flores2004; Hutchins & Marusic Reference Hutchins and Marusic2007; Reynolds et al. Reference Reynolds, Hayden, Castro and Robins2007). In rough-walled boundary layers, enhanced mixing can distort the instantaneous streaks, but they have still been observed in temporal- and composite-averaged fields (Coceal et al. Reference Coceal, Dobre, Thomas and Belcher2007; Barros & Christensen Reference Barros and Christensen2014; Fang & Porté-Agel Reference Fang and Porté-Agel2015). Persistent mean streamwise low-momentum streaks (manifesting themselves as weak streamwise vortices) have also been found in fully-developed canopy flows (Brunet et al. Reference Brunet, Finnigan and Raupach1994; Finnigan et al. Reference Finnigan, Shaw and Patton2009; Bailey & Stoll Reference Bailey and Stoll2013; Perret & Ruiz Reference Perret, Ruiz, Venditti, Best, Church and Hardy2013) (this will be further explored in § 5.4).

In the case of the mixing layer, these streaks are thought to originate from a secondary-instability mechanism which eventually causes vortices to become strained in the streamwise direction (Lin & Corcos Reference Lin and Corcos1984; Rogers & Moser Reference Rogers and Moser1992). When perturbations of streamwise vortices are in phase between adjacent rollers, they induce upward transport, resulting in streaks of low-momentum fluid. This in-phase formation of streamwise structures is commonly referred to as the ‘translative instability’ (Pierrehumbert & Widnall Reference Pierrehumbert and Widnall1982), which is depicted in figure 5 b). As stressed by Rogers & Moser (Reference Rogers and Moser1992), this is not exclusively a core instability (as in figure 5 b) but rather an instability of the whole flow which may contain sub-roller scales. Out-of-phase displacements, termed ‘helical pairing’, were also explored by Pierrehumbert & Widnall (Reference Pierrehumbert and Widnall1982). As depicted in figure 5(c), this mode of instability does not result in streamwise streaks because transport is vertically symmetric along a streamwise transect. Both the translative and helical pairing instabilities have been commonly observed during early evolution of the flow. Which of these two that dominates early stages the flow is sensitive to initial perturbations (Comte, Lesieur & Lamballais Reference Comte, Lesieur and Lamballais1992; Collis et al. Reference Collis, Lele, Moser and Rogers1994). The helical pairing modes are not unstable to high-wavenumber disturbances, while the translative instability is more broadband in nature (Pierrehumbert & Widnall Reference Pierrehumbert and Widnall1982). Pierrehumbert & Widnall (Reference Pierrehumbert and Widnall1982) investigated the persistence of these two instabilities into later flow development. They found that, although the helical pairing instability could be more amplified initially, it tends to die off with further flow development as pairing occurs and smaller scales emerge. Conversely, the translative instability was found to persist into later flow development, as it is still highly unstable for high-frequency perturbations. There is experimental evidence to suggest that it is the translative instability that dominates the late evolution of mixing-layer flow. The first is the observed phenomenon of streamwise streaks (Browand & Troutt Reference Browand and Troutt1980, Reference Browand and Troutt1985; Lasheras & Choi Reference Lasheras and Choi1988). Secondly, several authors have found that features of later flow evolution quantitatively agreed with the translative instability (Huang & Ho Reference Huang and Ho1990; Rogers & Moser Reference Rogers and Moser1992; Comte, Silvestrini & Bégou Reference Comte, Silvestrini and Bégou1998). For example, Huang & Ho (Reference Huang and Ho1990) found that spanwise length scales agreed with predictions of Pierrehumbert & Widnall (Reference Pierrehumbert and Widnall1982) for the translative instability. Rogers & Moser (Reference Rogers and Moser1992) found that growth rates of unstable modes corresponded to that of the translative instability. To our knowledge, no previous work has reported the dominance of the helical pairing instability for late evolution of a mixing layer (i.e. after full transition to turbulence has occurred).

Figure 5.

Illustration of three-dimensional vortex instabilities in mixing-layer flows (flow is from left to right). Vortex structures start as predominantly spanwise (a), then may evolve according to the translative instability (b) or the helical pairing instability (c). Shaded (unshaded) areas in (b) indicate regions of low-speed (high-speed) streaks in streamwise velocity.

Similarly, in boundary-layer flows, Klebanoff, Tidstrom & Sargent (Reference Klebanoff, Tidstrom and Sargent1962) originally noted the presence of a secondary instability that could lead to in-phase displacements of primarily spanwise structures corresponding to harmonic modes, which is usually referred to as a K-Type instability. This mode of instability has been reported to result in ‘peak-valley splitting’, which is qualitatively similar to what we described above for the translative instability of mixing layers (Klebanoff et al. Reference Klebanoff, Tidstrom and Sargent1962; Saric & Thomas Reference Saric, Thomas and Tatsumi1984; Herbert Reference Herbert1988). Boundary layers also have a helical pairing analogue with out-of-phase displacements corresponding to sub-harmonic modes (Herbert Reference Herbert1988). As in mixing layers, the selection of in-phase and out-of-phase instability is sensitive to initial conditions. Most studies have examined secondary instability in terms of its influence on the transition to turbulence. From a theoretical perspective, it is unclear which secondary-instability modes dominate the fully-developed boundary layer, however the presence of persistent streamwise streaks suggests a dominance of in-phase (harmonic) secondary instability.

In mixing layers, the spanwise wavelength of these streaks is slightly less than the primary roller spacing (about $(2/3){\it\Lambda}_{x}$ ; Pierrehumbert & Widnall Reference Pierrehumbert and Widnall1982; Huang & Ho Reference Huang and Ho1990). In boundary layer and canopy flows, the wavelength of the streaks is of the order of the boundary-layer height ${\it\delta}$ (Hutchins & Marusic Reference Hutchins and Marusic2007; Finnigan et al. Reference Finnigan, Shaw and Patton2009). We find from the start-up flow that the wavelength of the streaks is roughly equal to $L_{z}$ or ${\it\delta}$ . The study of Ghisalberti (Reference Ghisalberti2010) regarding the structure of an aquatic canopy also revealed what looks like in-phase displacements of spanwise structures, with a spanwise wavelength of about the boundary-layer height. Since spanwise length scales follow the same scaling for canopy and boundary layers, this would suggest that outer-layer instability is modulating the mean spanwise structure in the canopy layer. We find that the tallest hairpins tend to form in streamwise streaks of low streamwise velocity. Since the spanwise spacing between the streaks is consistent with instability in the overlying boundary layer, it would appear the secondary instability in the boundary layer creates a streaky structure that favours the formation of head-up hairpins. The work of Comte et al. (Reference Comte, Lesieur and Lamballais1992) also discussed how hairpin formation would be favoured in troughs of vorticity and inhibited in crests. This is consistent with observations in the ISL of hairpin packets coinciding with very long streaks in low-momentum fluid (Hommema & Adrian Reference Hommema and Adrian2003; Hutchins & Marusic Reference Hutchins and Marusic2007). As was also suggested by Hutchins & Marusic (Reference Hutchins and Marusic2007), it is conceivable that outer scaling modulates structure cycles in the roughness layer. We claim that alternating outer-layer patterns in the momentum field favours the evolution of predominantly spanwise structures near the canopy into hairpin structures that propagate into the overlying ISL. This process also agrees with Townsend’s attached-eddy hypothesis which states that inner and outer-layer structures are not independent.

5.1. Conditional averaging methodology

In order to educe the structure of a characteristic eddy, we largely followed the conditional averaging methodology of Finnigan et al. (Reference Finnigan, Shaw and Patton2009), which can be summarized as follows. We define a generic Eulerian variable ${\it\phi}(x_{i};t)$ as a trigger for conditionally sampling of the flow. The trigger variable was chosen such that when it is between some threshold range ${\it\alpha}_{min}$ and ${\it\alpha}_{max}$ , it presumably coincides with the presence of a turbulent structure of interest. To constrain the length scale of the structure, the trigger was only examined at a fixed vertical level of the flow $z_{T}$ . Thus at various instances of the flow, ${\it\phi}$ was searched for instances where ${\it\alpha}_{min}<{\it\phi}(x,y,z_{T};t)<{\it\alpha}_{max}$ . Local maxima of ${\it\alpha}_{min}<{\it\phi}<{\it\alpha}_{max}$ were identified, and their planar coordinates $\{x_{T},\,y_{T}\}$ noted. To help ensure threshold crossings were independent, smaller maxima within some specified radius of a larger local maxima were discarded. At each remaining local maximum, a realization of the flow was extracted ( $\tilde{u} ,\,\tilde{v},\,\tilde{w},\,\tilde{p},\,{\it\lambda}_{2}$ ), which was centred about $\{x_{T},\,y_{T}\}$ . Each realization was averaged to form a composite view of the flow when the threshold was exceeded.

In a similar manner as Finnigan et al. (Reference Finnigan, Shaw and Patton2009), we visualized the time evolution of the composite eddy structure. At each instance of a triggering event (time $t_{T}$ ), we also looked forward and backward in time at a fixed volume (centred in space at $\{x_{T},\,y_{T}\}$ and in time at $t_{T}$ ). The volumes were then composite averaged at times before and after $t_{T}$ to give the temporal evolution of the composite eddy.

5.2. Pressure maxima trigger

We begin by re-examining the composite averaging trigger of Finnigan et al. (Reference Finnigan, Shaw and Patton2009), which takes ${\it\phi}=\tilde{p}^{\prime }$ , with ${\it\alpha}_{min}=c_{min}u_{\ast }^{2}$ and ${\it\alpha}_{max}=c_{max}u_{\ast }^{2}$ ( $c_{min}$ and $c_{max}$ are positive constants). One reason this trigger was suggested by Finnigan et al. (Reference Finnigan, Shaw and Patton2009) is because previous work found that positive peaks in pressure generally coincide with scalar microfronts (Fitzmaurice et al. Reference Fitzmaurice, Shaw, Paw and Patton2004), which are presumably indicators of high flow coherence or structure.

Our methodology here differs slightly from Finnigan et al. (Reference Finnigan, Shaw and Patton2009) in that we visualize the total velocity field $\tilde{u} _{i}$ rather than the perturbation velocity field $\tilde{u} _{i}^{\prime }$ . Finnigan et al. (Reference Finnigan, Shaw and Patton2009) visualized their flow using perturbation velocity vectors, and ${\it\lambda}_{2}$ calculated from the perturbation velocity field ( $\tilde{u} ^{\prime },\,\tilde{v}^{\prime },\,\tilde{w}^{\prime }$ ). By using perturbation quantities to calculate ${\it\lambda}_{2}$ (we denote this as ${\it\lambda}_{2}^{\prime }$ ), Finnigan et al. (Reference Finnigan, Shaw and Patton2009) removed the mean vorticity $\partial \langle u\rangle /\partial z$ , and thus their composite structure corresponds to the perturbation vorticity field. Their approach has several potential pitfalls. One is that it could possibly remove a structure in which the mean vorticity plays a major role. Secondly, it can introduce spurious structures that are not present in the actual velocity field. We have chosen to use the total vorticity field in this section to eliminate these possible issues. However, we further explore the consequences of visualizing structures using the perturbation velocity and vorticity field in the Appendix.

Figure 8 gives the time evolution of our composite structure based on a positive pressure perturbation (and mean shear included in ${\it\lambda}_{2}$ ). The structure was obtained by following Finnigan et al. (Reference Finnigan, Shaw and Patton2009) and setting $c_{min}$ and $c_{max}$ respectively to 4.34 and 6.51, and $z_{T}/h=1.0$ . As noted by Finnigan et al. (Reference Finnigan, Shaw and Patton2009), varying the magnitude of $c$ within a reasonable range does not change the general shape of the structure. Rather we found it acted to add noise/distortion to the structure caused by not enough realizations ( $c_{max}$ too high) or including too many weak structures ( $c_{min}$ too low). An animation of the full time evolution of the composite structure can be found in supplementary movie 2.

Figure 8.

Composite eddy structure obtained by conditional sampling ( $\mathscr{N}=40,592$ realizations) based on spikes in positive pressure perturbation at $z/h=1$ . The structure is visualized using isosurfaces of ${\it\lambda}_{2}$ at a value of $0.01\cdot \text{min}({\it\lambda}_{2})$ . The isosurfaces are coloured according to fluctuations in the filtered vertical velocity. (a) $tu^{\ast }/h=-0.176$ , (b)  $tu^{\ast }/h=0$ , (c)  $tu^{\ast }/h=0.176$ , (d)  $tu^{\ast }/h=0.351$ . An animation of the full time evolution can be found in supplementary movie 2.

The composite structure emerges as a pair of predominantly spanwise structures located near the canopy top. The upwind structure primarily creates downward transport, while the downwind structure creates upward transport. The structure extends to about $\pm h$ in the spanwise direction before noise begins to break them up. As mentioned previously, we should not expect the composite rollers to extend the full width of the domain since irregularities will tend to disrupt their coherence far away from the trigger. This is especially true considering that the trigger location is between two structures. As time evolves, the structures remain at approximately the same spacing and shape and simply advect along the canopy top. At the time of the trigger, we find that the trigger location $\{x_{T},y_{T}\}$ is between the two structures, which is consistent with what was observed in the start-up flow.

Visual inspection indicates that the structures have a spacing of somewhere around $3h$ in the streamwise direction. Following Raupach et al. (Reference Raupach, Finnigan and Brunet1996), we estimated the theoretical mean streamwise structure spacing as ${\it\Lambda}_{x}=8.1L_{s}$ , where $L_{s}$ is the canopy-top shear length scale. Assuming the velocity profile approaches zero in the lower canopy, $L_{s}$ can be estimated as (Raupach et al. Reference Raupach, Finnigan and Brunet1996)

Using the profile in figure 7(a), we obtained a value of $L_{s}\approx 0.4h$ which gives ${\it\Lambda}_{x}\approx 3.25h$ . Thus, the composite structure roughly agrees with the expected structure spacing. This theoretical spacing will be further explored in § 5.2.1.

Our composite structure never has a strong three-dimensional shape for the duration of the composite evolution period that was examined. We can confirm that if we visualize the structure using the perturbation velocity/vorticity field (see figure 14 a), we indeed obtain a more three-dimensional structure that is similar to that reported by Finnigan et al. (Reference Finnigan, Shaw and Patton2009). In fact, this three-dimensional structure does not conflict with our conceptual model, and we do not contend that it is not a real feature of the flow. We find that the perturbation structure forms a ring vortex for most ${\it\lambda}_{2}^{\prime }$ isovalue choices. Finnigan et al. (Reference Finnigan, Shaw and Patton2009) also noted that certain isovalue choices showed a ring vortex. This perturbation structure agrees with the perturbation structure of a near-wall head-up hairpin shown by Kim & Moin (Reference Kim and Moin1986). It would appear that the three-dimensional structure is likely associated with a head-up hairpin as claimed by Finnigan et al. (Reference Finnigan, Shaw and Patton2009), but that removing the mean vorticity diminishes the surrounding two-dimensional structure.

We conclude that the three-dimensionality is largely associated with the perturbation vorticity field, but that a two-dimensional roller structure appears to dominate the pressure field and presumably the associated scalar microfronts/ramps. Again, this is not to say that the flow lacks substantial three-dimensionality. The composite structures themselves are a composite of many smaller three-dimensional vortices, which almost certainly undergo continuous localized interaction or pairing processes. The roller structure is also expected to have irregularities, including sinusoidal perturbations with wavelength of order ${\it\delta}$ resulting from mean streamwise streaks. We argue that it is an underlying quasi two-dimensional structure dominating the creation of observed pressure spikes and microfronts.

5.3. Reynolds-stress trigger

In this section, we created a composite averaging methodology to educe the composite structures associated with the generation of large Reynolds stresses. We recall from the start-up flow that the largest stress-generating structures (roller cores) coincided with regions of negative rather than positive pressure. It is reasonable then to consider a conditional averaging trigger based on pressure minima (cf. Hunt et al. Reference Hunt, Wray and Moin1988; Jeong & Hussein Reference Jeong and Hussein1995). However, in the fully developed very high Reynolds number case, the rollers are not traditional vortex tubes but rather a conglomeration of smaller vortices exhibiting rotational motion. Thus, we should expect the rollers to be associated with negative pressure on average, but not necessarily instantaneous pressure minima.

To directly seek out the dominant structures contributing to Reynolds stresses, we formulated a conditional averaging methodology based on a trigger of local maxima in $-u^{\prime }w^{\prime }=-(\tilde{u} ^{\prime }\tilde{w}^{\prime }+{\it\tau}_{13})$ . As also noted by Ganapathisubramani et al. (Reference Ganapathisubramani, Longmire and Marusic2003), one advantage of this trigger is that it does not bias toward either sweeps or ejections, since both can contribute to $-u^{\prime }w^{\prime }$ . The methodology we used is the same as described in § 5.2 except that a composite-averaging event is triggered based on local minima of regions where $c_{min}u_{\ast }^{2}(z_{T})<u^{\prime }w^{\prime }<c_{max}u_{\ast }^{2}(z_{T})$ , where here $u_{\ast }^{2}(z_{T})$ is (minus) the mean Reynolds stress at the trigger height rather than the canopy height. Threshold values of $c_{min}$ and $c_{max}$ were respectively given values of $-18$ and $-14$ .

Figure 10.

Time evolution of composite eddy structure obtained by conditional sampling ( $\mathscr{N}=25,345$ realizations) based on local minima of $u^{\prime }w^{\prime }$ at $z_{T}/h=1.5$ . The structure is visualized using isosurfaces of ${\it\lambda}_{2}$ at a value of $0.01\cdot \text{min}({\it\lambda}_{2})$ . The isosurfaces are coloured according to fluctuations in the filtered vertical velocity. (a)  $tu^{\ast }/h=-0.527$ , (b)  $tu^{\ast }/h=-0.351$ , (c)  $tu^{\ast }/h=0$ , (d)  $tu^{\ast }/h=0.439$ . An animation of the full time evolution can be found in supplementary movie 3.

Figure 11.

Composite eddy structure obtained by conditional sampling based on local minima of $u^{\prime }w^{\prime }$ at varying trigger heights of (a) $z_{T}/h=0.75$ , (b) $z_{T}/h=1.5$ , (c) $z_{T}/h=3.0$ . The structure is visualized using isosurfaces of ${\it\lambda}_{2}$ at a value of $0.01\cdot \text{min}({\it\lambda}_{2})$ . The isosurfaces are coloured according to fluctuations in the filtered vertical velocity.

The time evolution of the composite structure is shown in figure 10. An animation of the full time evolution can be found in supplementary movie 3. The composite structure emerges as a predominantly spanwise vortex structure that is at least $4h$ wide in the $y$ direction. The structure has the expected sense of rotation for a spanwise roller structure, as indicated by patters in $\tilde{w}^{\prime }$ . This predominantly spanwise structure supports the conceptual model Stages (a)–(b), and our previous hypothesis that the flow is dominated by an underlying two-dimensional structure. It is interesting to note that whether ${\it\lambda}_{2}$ or ${\it\lambda}_{2}^{\prime }$ is used to visualize the structure (see figure 14 b), a two-dimensional structure is still present. This is likely due to the fact that when the stress trigger is used, the structure is located at the trigger location. We can also confirm that if perturbation velocity vectors are used to visualize the structure (not shown), a two-dimensional structure is clearly present that is at least $4h$ wide in the spanwise. It is likely that the closer the structure of interest is located to the trigger location, the better defined it will be. The further the structure from the trigger location, the more it is liable to be eliminated by removing mean vorticity, and by composite averaging in general.

As time progresses, a hairpin-like structure forms above the spanwise structure inclined at approximately $45^{\circ }$ . This supports conceptual model stages (c)–(d). Away from the hairpin, the still predominantly spanwise portion of the structure generates sweep events, particularly on its downwind side. We do not observe any clear formation of a head-down hairpin structure, although one could argue that a ‘bulge’ exists where head-down structures would be located. Note also that the base of the structure is centred about $x=0$ at $t=0$ , in contrast to the structure in figure 8 where there was no structure in the immediate vicinity of $x=0$ at $t=0$ . This means that the actual structure is centred at the trigger location rather than surrounding it. As time progresses further, the head of the upper hairpin-like structure separates from the predominantly spanwise structure at the canopy top, supporting the ‘recycling’ stage of the conceptual model. Overall, the evolution of the composite structure is similar to what was observed in the start-up flow, and supports the conceptual model presented in § 4. It is impossible to fully support conceptual model Stage (b) with the composite structures, since their nature is to aggregate smaller structures. However, vortex breakdown is inevitable, and it is entirely inconceivable to imagine that a well-defined spanwise vortex tube exists given the extremely high Reynolds number and the frequently reported existence of an approximately $k^{-5/3}$ kinetic energy spectrum near the canopy (Finnigan Reference Finnigan2000; Bailey & Stoll Reference Bailey and Stoll2013).

The selected height of the trigger $z_{T}$ had an influence on the resulting composite structure and suggested that different structures dominate Reynolds stresses at different heights. Figure 11 shows the variation of the composite structure with $z_{T}/h=0.75$ , 1.5, and 3.0. At a trigger height of $0.75h$ , the composite structure is predominantly spanwise, and resembles the composite structure in figure 10 at times well before and after the triggering event. The structure has a slight bend around the trigger location, but as will be shown in the next section, this is because it likely resides in a low-speed streak. For a trigger height just above the canopy, we still obtain the predominantly spanwise structure, but an additional structure resembling that of a head-up hairpin emerges. As the trigger height is increased to several canopy heights, the resulting structure resembles the head of a head-up hairpin, although this is speculative since the details of the legs are not visible. This could be a deficiency of the composite-averaging method, which tends to lose structure detail far away from the trigger location. Qualitatively, this trend with varying $z_{T}$ agrees with the conceptual model presented in § 4.

5.4. Location of trigger events

If we examine temporally-averaged flow fields, we find that as in the start-up flow (figure 4 d) a streaked pattern emerges in many flow variables. Figure 12 shows temporally-averaged fields of streamwise velocity fluctuation $\overline{u^{\prime }}$ (figure 12 a), pressure fluctuation $\overline{p^{\prime }}$ (figure 12 b) and scalar fluctuation $\overline{{\it\theta}^{\prime }}$ (figure 12 c) at the canopy top. $\overline{u^{\prime }}$ has pronounced sinusoidal streaks in the streamwise direction, corresponding to weak streamwise vortices.

Figure 12.

Temporally averaged $x{-}y$ slices at $z=h$ of (a) streamwise velocity perturbation, (b) pressure perturbation, and (c) scalar perturbation. Spatial coordinates are normalized by the boundary-layer height ${\it\delta}=L_{z}$ . The right column shows the probability density function of conditional averaging trigger event spanwise locations for pressure trigger (squares) and Reynolds-stress trigger (circles). Note that the same plot is repeated for comparison with each $x$ – $y$ slice.

We find that the spanwise wavelength of our mean streamwise structures scales with $L_{z}$ , or equivalently the boundary-layer height ${\it\delta}$ . This would suggest that secondary instability in the overlying boundary results in in-phase displacements of primarily spanwise structures even near the canopy. The patterns in $\overline{u^{\prime }}$ and $\overline{{\it\theta}^{\prime }}$ suggest that certain streamwise regions of the flow favour upward transport.

To investigate the relationship with our detected coherent structures and mean streamwise structures, we created a probability density function (PDF) of trigger locations. This PDF is plotted alongside the temporally-averaged fields in figure 12(a–c). We found that most of the composite trigger events corresponded to mean streaks of negative $\overline{u^{\prime }}$ , positive $\overline{p^{\prime }}$ and positive $\overline{{\it\theta}^{\prime }}$ . Intuitively, it makes sense that since our composite-averaging methodology seeks out large Reynolds-stress events, most trigger events would happen in regions of low-momentum fluid which presumably favour the generation of head-up hairpin vortex packets. Not surprisingly, the resulting composite structure was associated with a head-up hairpin structure. We suspect that any conditional-sampling methodology seeking to capture large variance or flux events will be liable to bias toward these regions of low-momentum fluid and corresponding packets of hairpin vortices. Hence, they are likely to return a structure with some degree of three-dimensionality. This is exemplified by our Reynolds-stress trigger, and previously used methodologies based on empirical orthogonal function (EOF) decompositions (e.g. Finnigan & Shaw Reference Finnigan and Shaw2000). Care should be taken when using perturbation fields to visualize composite-averaged fields based on a point of large pressure, stress, variance, etc., because the perturbation field will tend to create a three-dimensional perturbation structure that surrounds the trigger point.

This result agrees with the start-up flow where we found that the tallest head-up hairpin structures tended to form along streamwise transects associated with low-momentum fluid. This also offers evidence to support our conceptual model, where we claimed that head-up hairpin formation is most commonly in phase between adjacent rollers.