The Role of Ambient Turbulence in Canopy Wave Generation by Kelvin–Helmholtz Instability

We test the hypothesis that internal waves observed in flow over forest canopies are generated by Kelvin–Helmholtz instability. The waves were observed at night, under stably stratified and weak wind conditions, with a horizontally scanning aerosol lidar and an instrumented tower. The lidar images are used to determine the salient wavelength and phase propagation velocity of each episode. Time series data measured at the tower are then used to form vertical profiles of background velocity and buoyancy just before each observed wave event. The profiles are input to the Taylor–Goldstein equation to predict the phase velocity, wavelength and period of the fastest-growing linear instability, and the results compared with the lidar observations. The observed wavelengths tend to be longer than predicted by the Taylor–Goldstein theory, typically by a factor of two. That discrepancy is removed when the theory is extended to account for the effects of ambient, small-scale turbulence.


Introduction
Mean flows in the atmospheric boundary layer are often non-inflectional, i.e. there is no local maximum in the shear except at the boundary, and inflectional instability is therefore not expected (Smyth and Carpenter 2019), an example being the log layer (Kaimal and Finnigan 1994;Tong and Ding 2020). But a forest or other vegetation layer reduces the shear near the ground, producing a shear maximum near the treetop height. Turbulent coherent structures often observed in forest canopies are therefore thought to result from inflectional, or more specifically Kelvin-Helmholtz, instability (hereafter KHI;Kelvin 1871;von Helmholtz 1890;Raupach et al. 1996;Smyth and Moum 2012). Here we test this hypothesis by comparing wavelengths, phase velocities and periods of a subset of 53 wave events observed over an orchard canopy during stably stratified conditions with those predicted via linear stability analysis of flows measured at a nearby meteorological tower. We find that the test is successful only when the theory is extended to account for the ambient turbulence that is present just prior to observed wave events.

Canopy Waves
Studies of turbulence in forest canopies over several decades have revealed coherent structures with length scales on the order of the canopy height (Raupach and Thom 1981;Finnigan 2000). In unstable or neutral conditions, these appear to begin as rollers and to rapidly evolve into arch-and hairpin-shaped structures (Bailey and Stoll 2016). In stably stratified conditions, that evolution is suppressed, and internal waves with parallel bands of crests and troughs perpendicular to the mean flow are observed (Mayor 2017, hereafter M17).
For example, dramatic oscillations with wavelengths of order 100 m were observed in time series data measured over 20-m tall forest canopies at Camp Borden, Canada (Lee et al. 1996) and in Prince Albert National Park, Saskatchewan, Canada . Regarding the mechanical origin of such oscillations, it has been hypothesized that "the turbulence and coherent motions near the top of a vegetation canopy are patterned on a plane mixing layer, because of instabilities associated with the characteristic strong inflection in the mean velocity profile." (Raupach et al. 1996). In other words, canopy waves are believed to be created via the mechanism of KHI, which we describe next.

Kelvin-Helmholtz Instability
We begin with some definitions. The local atmosphere is measured by a Cartesian coordinate system in which {x, y, z} are zonal (increasing eastward), meridional (increasing northward) and vertically upwards coordinates, {u, v, w} are the corresponding velocity components. The background horizontal velocity is denoted u h = {ū,v}. In our theoretical analyses, the overbar is taken to indicate a horizontal average. In the analysis of point data, as from a meteorological tower, it indicates a temporal average, here over a 5-min period.
Inflectional instability draws energy from the background flow at an inflection point of the background velocity profile or, more generally, at a local maximum of the shear magnitude: The instability mechanism is most easily understood as a constructive interference between vorticity waves that propagate horizontally above and below the inflection point and are therefore able to phase-lock (e.g., Carpenter et al. 2013;Smyth and Carpenter 2019;Liu et al. 2022). The resulting resonance creates exponentially growing plane wave structures in which departures from the mean state, e.g. u = u −ū, have the form: u (x, y, z, t) = û(z)e σ t+ι(kx+ y) whereû(z) is a complex vertical structure function, σ is the complex growth rate, subscript r indicates the real part, and ι = √ −1. The wave vector k = {k, } points horizontally. It may also be represented in polar form in terms of its magnitude κ = √ k 2 + 2 and its direction φ = tan −1 ( /k), measured counterclockwise from eastward. The wavelength is λ = 2π/κ, the phase propagation speed is c = −σ i /κ and the period is τ = λ/c.
To optimize the inflectional instability process, the wave vector aligns with the background shear {ū z ,v z } at the shear maximum. The phase speed lies between the minimum and the maximum ofũ(z), the component of the background velocity profile in the direction of the maximum shear defined above (the semicircle theorem; (Howard 1961)). The wavelength of these plane waves scales with the thickness of the shear maximum where the instability is focused. For example, in an idealized shear layerū ∝ tanh(2z/h), where h is the layer thickness, the wavelength is approximately 7h. The scale factor 7 can vary considerably, though, when the velocity profile has a more complicated form or when other effects are active.
KHI is the variant of the inflectional instability that occurs in a stably stratified fluid. That state is defined by the positivity of the squared buoyancy frequency, computed here as: where g is the gravitational acceleration (here taken to be 9.81 m s −2 ), T v is the background virtual temperature and v is the corresponding potential temperature. The gradient Richardson number is a dimensionless number that quantifies the strength of the shear relative to the stratification: The growth rate of the inflectional instability vanishes if the stratification is strong enough that the minimum gradient Richardson number, min z {Ri}, exceeds 1/4, i.e. the buoyancy frequency N > S/2 at all z (Miles 1961;Howard 1961). On the other hand, if the stratification is weak enough, the instability grows. The wave vector is aligned with the background shear vector as in unstratified flow, and the wavelength is only slightly modified. At large amplitude, the instability is manifested in trains of overturning billows that we call KH billows or (casually) KHI.
Stably stratified shear flows also support constructive interference involving gravity waves, leading to distinct mechanisms such as the Holmboe and Taylor-Caulfield instabilities (Smyth and Carpenter 2019;Arnqvist et al. 2016). These will not be considered here.
Because it grows on parallel shear flows with min z {Ri} < 1/4, we associate KHI with low Ri, i.e. strong shear and/or weak stratification. Moreover, as the instability grows to finite amplitude, it suffers a sequence of secondary instabilities that lead to a fully turbulent state Peltier 1985, 1991;Caulfield and Peltier 1994;Mashayek and Peltier 2012a, b). We therefore expect to observe enhanced turbulence associated with KHI. As we will see, the observed canopy waves behave differently: they appear at times when Ri is not especially small, and turbulence is weak.
Stratified shear flows are subject to both viscous and diffusive effects ). These need not be molecular in origin; they can also result from turbulence on scales much smaller than the instability, possibly left over from previous instabilities Kaminski and Smyth 2019). While the classical KHI model describes an infinitesimal disturbance growing on an otherwise parallel shear flow, natural flows almost always exhibit some degree of small-scale turbulence. Viscosity, either molecular or turbulent, tends to damp KHI in two ways: by directly retarding the unstable motions and by diffusing the background shear that drives them. Diffusion acts similarly to smooth out the instability, but its effect on a stably stratified background flow is to reduce the stratification, thus enhancing the growth of the instability. Ambient turbulence generally varies with height, which can lead to more complex effects (Li et al. 2015).

Are Canopy Waves Created by KH Instability?
Because the inflectional instability mechanism results from constructive interference between waves above and below a shear maximum, it cannot occur when the shear maximum is at the ground, as in a standard logarithmic wind profile. However, vegetation can mix the background shear near the ground, creating an inflection point at which instability may grow. This suggests a possible explanation for the observed canopy waves (Raupach et al. 1996).
Analysis of the Saskatchewan canopy wave observations indicated that the waves were generated by shear near the top of the canopy (Lee 1997). Moreover, the observed waves shared some identifying features of KHI: their phase speed and direction were equal to the background wind near the center of the shear layer, and their amplitude decayed vertically away from that layer. Linear stability analysis confirmed that the observed profiles of velocity and density supported unstable modes. The computed growth rate, however, was not a maximum at the observed wavelength as would be expected.
A 20-m forest in Sweden produced inflectional flows and canopy waves that were welldescribed by KHI theory, even though tree drag and heat transfer were neglected (Arnqvist et al. 2016). In particular, the fastest-growing mode was identified and found to match the observed frequency. The interpretation suggested Holmboe-type instabilities (Sect. 1.2).
Although the connection between canopy waves and KHI is theoretically plausible and is consistent in many respects with existing observations, other processes, such as streamline displacement via wave-turbulence interactions (Sun et al. 2015) and the Jeffrey drag mechanism (Pulido and Chimonas 2001), can create waves in the forest environment. These mechanisms are alternatives to KHI for the creation of canopy waves.

The Canopy Horizontal Array Turbulence Study
The Canopy Horizontal Array Turbulence Study (CHATS; Patton et al. 2011) was conducted in 2007 to observe the fine-scale structure of turbulent flow over an idealized forest of uniform height and density. The experiment was conducted in a 1.6 km 2 walnut orchard on flat terrain in Dixon, California, from mid-March through mid-June. The Raman-shifted Eye-safe Aerosol Lidar (REAL; Mayor et al. 2007) was deployed at a range of 1.61 km from a 30-m tall instrumented tower installed in the orchard. The average height of the trees was about 10 m and the nearly-horizontal lidar scans intersected the tower at an altitude of about 18 m above ground level.
Most observations of KHI in the lower atmosphere are made with radars (Richter 1969;Browning and Watkins 1970;Gossard et al. 1970;Chapman and Browning 1997), instrumented towers (Einaudi and Finnigan 1981;Finnigan et al. 1984;De Baas and Driedonks 1985), and sodars (Emmanuel et al. 1972;Hooke et al. 1973;Emmanuel 1973;Petenko et al. 2012). The observations are often limited to vertical profiles. To obtain 2D spatial data, which is necessary for the most reliable measurement of wavelength, scanning remote sensors are required. For example, Blumen et al. (2001) and Newsom and Banta (2003) describe an episode of turbulence associated with KHI observed in the stable boundary layer using a scanning Doppler lidar. In that work, the Doppler lidar scans provided vertical cross-sections of the waves from one event. In our work (M17), the REAL was able to observe the waves on near-horizontal crosssections and to capture 53 events over three months (Fig. 1). Graphical analysis of the images resulted in a wavelength assessment for 52 of these cases. The power and fast scan update rate of the REAL enabled the determination of the phase propagation velocity in 39 episodes where the lidar frame rate was fast enough to unambiguously determine wave crest displacement. We are not aware of any other atmospheric observational datasets that provide horizontal images of canopy wave spacing and phase velocity and coincident vertical profiles of thermodynamic and wind data. Herein we use these data to see how the prevailing KHI theory can be modified to link the observed vertical profiles and the wave characteristics revealed by the lidar.
During the 53 events identified by M17, a shear maximum was almost always found near z = 10 m, the nominal treetop height. Most of the horizontal lidar images intersect the tower (and therefore the vertical profiles) near 18 m height. Plane wave crests were nearly perpendicular to the background shear, consistent with the KHI theory outlined above. Moreover, waves propagated in the direction of the wind, and with about half the maximum wind velocity, consistent with the semicircle theorem (M17). Observed wavelengths were not evidently related to the length scale of the background shear, as is the case with KHI. As a first estimate, we may take the thickness h of the shear layer to be 3-5 m (M17), in which case the KHI wavelength estimate 7h becomes 21-35 m. In fact, based on the lidar images, the observed wavelengths ranged from 30 to 100 m, considerably longer than the theoretical estimate. That disparity is the primary motivation for the present study. We will argue that it results from ambient turbulence, which greatly affects the preferred wavelength of KHI.

Outline
The discussion is organized as follows. In Sect. 2 we composite the tower data to examine the evolution of the background flow, waves and turbulence over a diurnal cycle (Sect. 2.1) and over a canonical canopy wave event (Sect. 2.2). In Sect. 3 we use profiles of observed shear and density from the CHATS tower data, together with ambient turbulence having various assumed amplitudes, to calculate wave characteristics such as wavelength and phase velocity based on linear stability analysis. We then compare the results with observations from REAL. Conclusions are summarized in Sect. 4.

Evolution of the Mean Flow, Waves and Turbulence
In this section we examine the time-dependent relationship between KHI, canopy waves, turbulence and the evolving background state using composite representations of a diurnal cycle and of a canopy wave event. We also take an initial look at turbulent, or eddy, viscosity based on the tower data.

Fig. 3
Nocturnal probability density functions for the logarithm of a the minimum value of Ri and b the vertical velocity variance for 9 < z < 18m and 19:00-07:00 local time. c Values measured over 5-min intervals immediately preceding canopy wave observations

Flow Properties Prior to Canopy Wave Events
For each wave episode registered by REAL, a mean (or "background") state is defined by averaging tower measurements of wind velocity, temperature and humidity at over a 5-min period immediately prior to the appearance of the waves. In these states, the shear magnitude S(z) is a maximum near z = 10 m, coinciding with the nominal treetop height (Fig. 2a, blue curve). This height is also a local minimum of the squared buoyancy frequency N 2 and of Ri (Fig. 2b). The minimum of Ri is less than 1/4 just prior to all wave events (Fig. 3), consistent with KHI theory (Sect. 1.2). The layer in which Ri < 1/4 extends from slightly below z = 10m to approximately z = 16m.
The probability density function of Ri for nighttime data is unimodal with mode near 0.04 (Fig. 3a). Also unimodal was the vertical velocity variance w w (Fig. 3b). Canopy waves emerged only when Ri was less than 1/4 and w w was relatively small (Fig. 3c), indicating minimal small-scale turbulent transport despite the potential for KH instability.

The Composite Diurnal Cycle
Properties of the composite day (Figs. 4, 5, and 6) are educed by compositing the tower data recorded during each hour of the day for the 3-month duration of CHATS. For context, we begin by examining the upper level wind velocity, averaged over sensors at z = 23 and 29 m and over each instance of each hour. The meridional component (Fig. 4a, red curve) shows a diurnal variation with northerly flow from 04:00 to 14:30 LT and southerly flow at other times. This reversal coincides approximately with minima in the resultant averaged wind speed (blue). The zonal wind (Fig. 4a, yellow) is westerly and is relatively weak and steady.
The cause of the meridional wind reversal is the diurnal heating and cooling of California's Central Valley relative to the steadier marine airmass over the Pacific Ocean 100 km to the west (Zaremba and Carroll 1999). The California Coast Ranges separate the cool marine air from the diurnally changing inland air except for a gap in the San Francisco Bay area. Air entering the Central Valley through the gap (approximately 55 km to the southwest) moves northward over Dixon during the afternoon and evening until late night and early morning cooling in the Central Valley relaxes the horizontal pressure gradient enabling a northerly prevailing flow. The prevailing northerly flow is the result of the north-south orientation of the Central Valley and synoptic scale flow around a semi-permanent anticyclone that exists over this region.  Evaluating the vertical derivatives ∂ū/∂z and ∂v/∂z at the height where the shear magnitude is a maximum (usually near 10 m), and repeating the same compositing operations described above, we have the diurnal cycle of the shear magnitude and the zonal and meridional shear components (Fig. 4b). These vary similarly to the upper level winds: ∂ū/∂z (yellow curve) is weak and approximately constant while the dominant component ∂ū/∂z (red) changes sign in the early morning and the early evening. The resultant shear magnitude (blue) decreased monotonically from midnight to sunrise, reaching an early morning minimum roughly coincident with the meridional shear reversal.
The squared buoyancy frequency N 2 (Fig. 5a), evaluated at the shear maximum, is considerably smaller than S 2 /4, suggesting a strong likelihood of instability (Miles 1961;Howard 1961). During the day, N 2 dropped to near or below zero, indicating statically neutral or unstable conditions. Underlying the S 2 and N 2 curves on Fig. 5 is a histogram of the number of wave events observed in each hour of the REAL imagery. All observed events occurred at night, with event frequency increasing to a maximum at 06:00, just before sunrise.
This composite diurnal cycle presents some surprises. Ordinarily, we associate KHI with strong shear and weak stratification, or low Ri. But the present observations suggest a different scenario, in which instability is detected when shear is relatively weak, stratification is strong and thus Ri is higher than usual. Reflecting the variation of the shear, Ri is low during the day, rises rapidly around sunset and then more gradually to a maximum before dropping rapidly around sunrise (Fig. 5b, blue curve).
Contrary to the usual KHI model, wave events are not associated with the lowest Ri. Instead, the Miles-Howard criterion min z {Ri} < 1/4 is virtually always satisfied, usually by a large margin. In this near-homogeneous state, shear instabilities evolve rapidly and chaotically, transferring energy both downscale (as turbulence) and upscale (via vortex pairing; Winant and Browand 1974;Klaassen and Peltier 1989). Any recognizable monochromatic wave signal is quickly overwhelmed. Only when Ri increases to a moderate value can waves survive long enough to be detectable. Thus, waves are always detected at night and are more common toward sunrise when Ri is greatest.
The turbulent heat flux (Fig. 6a) was positive (upward) from approximately 07:00 to 18:00 LT and negative at other times. 1 The vertical velocity variance (Fig. 6b) increased steadily throughout the day, decreased rapidly around sunset, then decreased more gradually through the night. The peak event frequency at 06:00 coincides with the minimum vertical velocity variance, i.e., the least-turbulent ambient conditions. These observations suggest that the vertical motions revealed in the velocity variance and the heat flux ( Fig. 6) are dominated not by the wavelike instabilities detected by REAL but rather by turbulence.

The Composite Canopy Wave Event
The composite event is formed by shifting the time during each canopy wave observation so that t = 0 when the waves are first visible in the REAL imagery. We then average over all 53 episodes and show the results from one hour before to one hour after the time of first visibility.
The upper level wind (Fig. 7a) does not change much in strength prior to wave emergence, though it changes in direction. The maximum shear (Fig. 7b, black curve, 7c, blue) trends downward, showing only a slight increase just before the event is observed. During the interval from − 30 to − 5 min., S 2 decreases by 31%. The buoyancy frequency changes very little (red); the increase from − 30 to − 5 min. is 9%. Decreasing shear and increasing stratification combine to create a 37% increase in Ri. In summary, the variables that normally govern inviscid KHI do not indicate any systematic change that would account for the emergence of instability at the particular time when it is observed. The change in the median long-wave radiative heat flux (Fig. 7e) is only ∼5%. Together with the relative steadiness of N 2 and Ri, this indicates that destabilization by radiative heating (Cava et al. 2004) is not a critical factor in these observations.
Of the variables examined here, the vertical velocity variance (Fig. 7f) shows the only change that could explain the appearance of instability, decaying by 32% prior to wave detection. We interpret this to suggest that, prior to wave emergence, instability is suppressed by ambient turbulence, and appears only when the latter has decayed sufficiently to permit sustained growth. The variance then increases rapidly as the wave signal grows.

Observational Estimates of Eddy Viscosity and Diffusivity Prior to Canopy Wave Events
The classical KHI theory is likely to be insufficient to explain these observations because the flow is turbulent before the instability sets in. The effect of this pre-existing turbulence is modeled most simply by choosing values for the corresponding eddy viscosity and turbulent diffusivity just prior to the onset of KHI. Vertical eddy viscosity A v (z) and diffusivity K v (z) were calculated from the tower data via the relations: where primes indicate departures from the 5-min mean. During the 5-min period preceding each of the 53 observed events, these values are typically between zero and 0.2 m 2 s −1 (Fig. 8a). Above the treetops (z ≈ 10 m), median values of eddy viscosity and diffusivity are similar, i.e. the turbulent Prandtl number is near unity, as is commonly observed in strong turbulence (the Reynolds hypothesis). Both coefficients decay smoothly upward from z ≈10 m. Below z ≈ 9 m, A v and K v are unequal and highly variable. It is likely that, in the well-mixed air near the ground, A v and K v are ratios of small quantities and are therefore Is this ambient turbulence sufficient to affect KHI? We will determine this explicitly later, but for now we consider a simple estimate based on the assumption of uniform A v . Growth is affected significantly when the Reynolds number Re = h u/A v is smaller than ∼ O(10 2 ), and is arrested entirely when Re < O(10) (Maslowe 1973;Smyth et al. 2013). If we take a typical velocity scale u to be 1m s −1 and the layer thickness to be 4 m, then Re = 100 corresponds to A v = 0.04m 2 s −1 , similar in magnitude to the median measured value at treetop height quoted above. The condition Re < 10, required to stabilized KH instability, requires A v > 0.4m 2 s −1 . We therefore expect that eddy viscosity associated with ambient turbulence (Fig. 8) will have a significant damping effect on the growth of KH instability, but will not stabilize it entirely.
In the remainder of the paper, we demonstrate that ambient turbulence has a significant effect on the growth of the observed instabilities. For simplicity in this initial investigation, we assume that ambient turbulence is stationary, homogeneous and isotropic.

Comparison of Observed Canopy Waves with KHI Via Linear Stability Analysis
We now ask whether the observed canopy waves could have been generated by KHI. To answer this, we perform explicit linear stability analysis of the mean profiles measured just before each wave event. We predict the phase speed, direction and wavelength for each event and compare the results with the available observations from REAL. We also explore the possibilities that KHI is modulated by pre-existing turbulence or by air-vegetation fluxes.
The following two subsections give details of our methods for deriving wave parameters via the REAL imagery (subsect. 3.1) and via linear stability analysis (subsect. 3.2). Results of the comparisons are in subsect. 3.3.

Estimating Wave Parameters from REAL Imagery
In Fig. 1 we show 1 km 2 areas from horizontal scans of the REAL from four sample canopy wave episodes. A full description of how the lidar data were collected and processed and the wave characteristics were extracted is contained in M17. To summarize the final steps, wave characteristics were obtained graphically by projecting the images on to a whiteboard and hand tracing the wavecrests. Wavelength was determined by measuring the average distance between two crests in the vicinity of the tower. Phase propagation speed was determined by observing the displacement of a wave crest over sequential scans, while the direction was by definition normal to the wave crests. Figure 1 reveals the significant variability that exists among cases, as well as the spatial variability that may occur in a single image. In some cases, if the waves were not clear at the location of the tower, the measurement was made as close to the tower as possible. If the wave characteristics changed over time, an average was chosen. The goal of the analysis was to assign a single value of wavelength and phase velocity for each episode. A first attempt to determine a single wavelength and orientation per image objectively by computer algorithm was recently reported by Mifsud et al. (2021) and the results support the subjective values used in M17 and herein.

Linear Stability Analysis
KHI is obtained as a solution of the Taylor-Goldstein (hereafter TG) equation, which describes normal mode disturbances (defined in Sect. 1.2) evolving on an inviscid, stratified, parallel shear flow. The theory has been extended to include viscosity and diffusivity due to ambient small-scale turbulence (Liu et al. 2012;Lian et al. 2020), and to approximate the momentum and heat fluxes between the flow perturbation and the vegetation (Lee 1997). The Appendix describes the basic theory and methodology; here, we describe modifications unique to the present project.
For small perturbations from the background flow, the momentum exchange with the trees (e.g.,Lee 1997) is approximated by an acceleration −C d L f √ū 2 +v 2 u . The constant C d is a drag coefficient, typically set to 0.15, while L f (z) is the leaf area density. Based on the CHATS measurements (Patton et al. 2011), we assume a piecewise-linear form for L f : with tree height h = 10 m and z max = 8m. The rate of heat transfer from the trees is modeled similarly as −C h L f √ū 2 +v 2 b , where b = −gρ /ρ is the buoyancy perturbation and C h = 0.10.
We assume that the ambient turbulence is homogeneous, isotropic and has unit Prandtl number so that the viscosities and diffusivities that work on vertical and horizontal gradients, #2 Next, we add fluxes of momentum and heat between the air and the vegetation with coefficients C d = 0.15 and C h = 0.10 following Lee (1997). #3 Third, we retain C d = 0.15 and C h = 0.10 but add eddy viscosity/diffusivity with a fixed value A 0 = 0.06m 2 s −1 , a compromise between the observational estimates of A v and K v at treetop height (Fig. 8). #4 Finally, we determine for each case the value of A 0 that maximizes agreement between the computed wave vector of the fastest-growing mode and the estimate based on the REAL imagery (see subsect. 3.1). We repeat the calculation of the fastest-growing mode for a range of values of A 0 between zero and 0.4, then choose the A 0 that minimizes the difference in wave vectors. This model also retains C d = 0.15 and C h = 0.10.
We solve these models using as input profiles of horizontal velocity and buoyancy obtained from the CHATS tower measurements. For each of the 39 observed events for which wave parameters have been measured, the corresponding mean state is identified as the velocity and buoyancy profiles obtained by averaging in time over the 5-min period immediately preceding the first appearance of the billows in the REAL imagery. The 5-min averaging interval is intended to be long enough to average out the slowest observed wave periods (60 s) yet short enough to preserve fluctuations in the mean flow that create instability. To maintain the validity of the normal mode solution form, only modes with growth rate faster than 1/(5 min)=0.033 s −1 are retained. Wind above the uppermost measurement at 29 m is assumed to vary linearly, thus ensuring that U zz = V zz = 0 and avoiding spurious shear instabilities. The buoyancy frequency is assumed to be constant above 29 m.
Generalizing the approach of Lee (1997), we make no assumption about the wavenumber and propagation direction but instead scan over those parameters to identify the fastestgrowing instability. Wave properties, including speed and direction of phase propagation, wavelength and period are compared with estimates based on the REAL observations described in M17 (details in subsect. 3.1). By investigating not just a single case but instead 39 cases in which waves were observed, we are able to draw general conclusions about the validity of the linear model.

Results
We first describe a representative case that illustrates the methods and the main results. We then generalize the results via a statistical analysis of 39 events for which clear observations of canopy wave parameters are available.

Case Study: Event #40
To illustrate the analysis, we first examine event #40 (in the numbering of M17), which began at 03:33:34 LT (local time is URL -7 h) on 2 June 2007 and lasted for 6 min, 10 s. The background velocity profile (Fig. 9a) may be crudely described as a shear layer between ∼ 5 m and ∼ 23 m in which the wind increased from near zero to ∼ 2 m/s from the WSW.
The mean squared shear profile (Fig. 9b, blue curve) shows a sharp maximum at 10 m. The squared buoyancy frequency (Fig. 9b, red curve) is generally positive but is, importantly, less than S 2 /4 (implying Ri < 1/4) between 8 m and 23 m. The shear maximum at 10 m is accompanied by a stratification minimum (suggesting previous mixing) and thus very low Ri. The flow is therefore potentially unstable to shear-driven instabilities.
We now solve the stability problem for a range of wave vectors {k, } and plot the growth rate σ of the fastest growing mode at each point (Fig. 9c). For this illustrative calculation, we c Growth rate versus zonal and meridional wave vector components. Blue arrows indicate the observed wave vector; black arrows highlight the wave vector of the fastest-growing instability, κ = 0.2 m −1 , φ =50 • . d The shear production rate (normalized by its maximum value) use the TG model (#1 in Sect. 3.2). Instability is found over a broad range of wave vectors, with maximum growth rate 0.033 s −1 at κ = 0.2 m −1 (i.e. wavelength 30 m) and direction 50 • north of east. At this growth rate, the disturbance grows by a factor e in about 30 s. The period in which the waves are observed to grow includes several of these e-folding intervals, showing that the growth rate is sufficient to explain the observations. Note also that the time scale for instability growth is much faster than the 5-min averaging interval that defines the background flow. We are therefore justified in using the normal mode solution form (2). If this were not true we would reject the result.
The lower half of Fig. 9c is a reflection of the upper half due to the symmetry of the equations. The black arrow indicates the fastest-growing mode over all {k, l}. The indeterminacy in the wave vector of the fastest-growing mode due to the aforementioned symmetry is resolved by requiring that the phase propagation speed c be positive. The shear production rate (Fig. 9d) for this mode is sharply peaked near z = 10m, confirming that the instability draws its energy from the shear just above the treetops.
The blue arrow in Fig. 9c indicates the wave vector as measured graphically from the REAL imagery. Comparing the arrows, we see that the linear analysis predicts the direction of the wave vector quite accurately, but overestimates its length. More specifically, the observed wavelength of the billows is about twice that predicted by the stability analysis.
The calculation shown in Fig. 9 has been simplified by the omission of ambient turbulence effects and drag and heat transfer due to the leaves (i.e., model #1). We now ask how the result will change when those effects are accounted for. We first add the vegetation terms as specified in model #2, with coefficients C d = 0.15 and C h = 0.10. The result is shown in Fig. 10a. The main difference is an overall reduction in the growth rate; the maximum drops from 0.033 s −1 in the TG limit to 0.027 s −1 with the transfer terms included. The wave vector of the fastest-growing mode is virtually unchanged.
We next add eddy viscosity and diffusivity, with A 0 = 0.03m 2 s −1 . The growth rate distribution now changes markedly (Fig. 10b). (Note that the axes have been altered from Fig. 9 to make the central region of the plot more visible.) At large wavenumbers (around the outside of the diagram), growth rates are greatly reduced. The two large regions of instability have Fig. 10 Growth rate versus zonal and meridional wave vector components for event #40 with eddy viscosity/diffusivity amplitude A 0 = 0 (a), 0.03 m 2 s −1 (b), 0.06 m 2 s −1 (c) and 0.17 m 2 s −1 (d). Arrows identify the fastest-growing unstable mode (black) and the observed wave vector (blue). In all cases C d = 0.15 and C h = 0.10 moved inward, indicating that smaller-wavenumber features are less vulnerable to viscous damping. The fastest-growing instability (black arrow) now has a wave vector quite similar to the observation (blue arrow).
We now repeat the experiment with A 0 increased to 0.06 m 2 s −1 (model #3), and then to 0.17 m 2 s −1 . As shown in Fig. 10c and d, the region of fastest growth continues to move toward smaller wavenumbers (i.e. longer waves). The final eddy viscosity shown, 0.17 m 2 s −1 , yields a growth rate of 0.03 s −1 , beyond which the time scale for instability growth is longer than the 5-min time scale of the background flow and the normal mode solution is no longer valid.
In Fig. 11, the wave vector predicted for A 0 = 0.03m 2 s −1 is shown by an arrow superimposed on the REAL image for event #40. The reader may verify the correspondences in both magnitude (or wavelength) and orientation. Figure 12 gives a more thorough depiction of the dependence of wave properties on the assumed eddy viscosity A 0 . The growth rate (Fig. 12a) decreases monotonically with increasing A 0 . The wavenumber κ (Fig. 12b) decreases sharply as A 0 is increased from zero, then becomes independent of A 0 at higher values. The dependence of φ on A 0 (Fig. 12c) is less consistent, but the prediction generally remains within a few degrees of the observed value. The phase speed (Fig. 12d) increases rapidly then flattens out at a value close to the observations. The discrepancy in the wave vector (Fig. 12e) shows a clear minimum at A 0 = 0.03m 2 s −1 , the value used for Fig. 10b. The discrepancy in the phase speed (Fig. 12f) also drops rapidly with the introduction of eddy viscosity, decreasing by an order of magnitude at A 0 = 0.03m 2 s −1 .
The stability boundary, at which the growth rate drops to zero due to strong eddy viscosity, is of interest although it is not formally valid as it violates the frozen flow assumption. In Fig. 13 Growth rates from the linear stability analysis of 39 observed wave events. Models including fluxes due to vegetation and ambient turbulence are compared with the original TG equation: a model #2, with air-tree fluxes following (Lee 1997), b model #3: air-tree fluxes plus turbulence with eddy viscosity/diffusivity A 0 =0.06 m 2 s −1 , c model #4: air-tree fluxes plus turbulence with A 0 fitted to maximize agreement with the observed wave vector. d Growth rate from model #4 versus fitted eddy viscosity/diffusivity A 0 . For the point marked , the optimal value A 0 = 0 does not fit on the logarithmic scale this case, the eddy viscosity needed to bring the growth rate to zero is A 0 = 0.26m 2 s −1 (Fig. 12a).
The difference in effect between tree fluxes and turbulence is not unexpected. The tree fluxes are modeled as Rayleigh damping terms, which have no scale dependence. In contrast, the turbulence effects are modeled using Laplacian operators, which act more strongly on small-scale structures.

Statistical Analysis of 39 Events
We now apply our four LSA models described in Sect. 3.2 to the full ensemble of 39 canopy wave events for which subjective estimates of the phase propagation speed, direction. wavelength and period are available (M17). The TG model gives growth rates between 0.018 s −1 and 0.083 s −1 , with a median value 0.036 s −1 . The addition of trees (model #2) reduces the  Table 1 gives the corresponding comparisons in quantitative form median growth rate to 0.032, a reduction of 12% (Fig. 13a). The further addition of viscosity with A 0 = 0.06m 2 s −1 (model #3; Fig. 13b) creates a general decrease in growth rate, with median reduced by more than a factor of two. Replacement of the fixed value A 0 = 0.06m 2 s −1 with the individually fitted values (model #4; Fig. 13c) reduces the median growth rate by a further factor of two. The optimal eddy viscosity values delivered by model #4 range from zero to 0.3 m 2 s −1 , with larger values tending to produce smaller growth rates (Fig. 13d).
The optimal values of A 0 have median 0.06 m 2 s −1 , and quartiles 0.03 and 0.14 m 2 s −1 . These ranges are strikingly similar to the range of measured values at the shear maximum (Sect. 2.4). However, there is little correlation between measured and inferred values for individual cases.
We next compare the wave vector, phase speed and period computed using the four linear models with the subjective estimates of M17. The results shown previously for case #40 (Sect. 3.3.1) turn out to be generally true. We begin by discussing the TG, with reference to the top row of Fig. 14, and the third column of Table 1. The direction of propagation of the fastest-growing instability is reasonably consistent with that of the observed wave (Fig. 14a). This is close to the direction of maximum shear, consistent with KHI theory (e.g., Smyth and Carpenter 2019). The median of the observed directions (Table 1, row 4) is 65 • north of east, whereas the median of the LSA predictions (row 5, column 3) is 85 • , i.e. the prediction is high by 20 • , or 31% relative to the observed value (row 6). By coincidence, the median discrepancy in individual cases (computed as median absolute difference between predicted and observed directions; row 7) is also 20 • .
The approximate agreement found for the propagation direction is not found in the case of the wavenumber κ (Fig. 14b, Table 1, column 3, rows 8-11). Instead, both the difference of medians and the spread are greater than a factor of two, as was illustrated previously for case #40 (Fig. 9c). Points clustered at the top of Fig. 14b have κ = 0.4m −1 , the largest Table 1 Comparison of observed and computed wave scales using four LSA models Differences are also given, in parentheses, as a percentage of the observed median value included in the analysis. The true value is even larger, i.e. further from the observed wavenumber.
The median phase propagation speed c (Table 1, rows 12-15) is low by 18% while individual cases typically differ by 40%. The period (rows 16-19) is systematically low by 48% while individual values differ by 58%.
To summarize: the TG model reproduces the direction well and the speed to within a few tens of percent, both in the median and in individual cases, but it overestimates the wavenumber (and underestimates the period) by more than a factor of two. This suggests that canopy waves are not well represented by KHI, at least not in its classical, inviscid form. We now ask whether the agreement can be improved by accounting for air-tree fluxes and ambient turbulence using the parameterizations described above.
Adding the tree fluxes alone (model #2) does not improve the correspondence with observations (Table 1, column 4). In particular, the extreme discrepancy in the wavenumber is unchanged. In other cases, the discrepancies are actually slightly greater.
We turn next to the effect of eddy viscosity/diffusivity due to ambient turbulence. To begin, we add turbulent viscosity and diffusivity at the level of A 0 =0.06 m 2 s −1 , approximately the median of the measured values at treetop height (Fig. 8), i.e. model #3. Results may be found in Fig. 14, third row, and column 5 of table 1. All discrepancies are now significantly reduced. In particular, the discrepancy in wavenumber is reduced from a factor of two to a few tens of percent (table 1, rows 10,11).
These promising results conceal a serious deficiency. In 6 of the 39 cases studied, the stability analysis failed to yield a growing mode. The highest growth rate was either less than the cutoff value 0.0033 s −1 or negative, indicating that the value A 0 =0.06 m 2 s −1 is beyond the stability boundary. This is a serious failure since the REAL images clearly show an instability. There are multiple possible reasons for this. For example, the waves may have grown slowly enough that their properties were determined by events more than 5 min previous to their becoming visible. Alternatively, the variation of eddy viscosity and diffusivity with height, neglected here, may be important. But the most likely explanation, is that the single value A 0 =0.06 m 2 s −1 is not applicable in all cases, and is even beyond the stability limit for some. We now seek to remedy that problem.
In model #4, we employ for each case the value of A 0 that optimizes the correspondence in the wave vector. When the stability analysis is conducted using these optimal values (table 1, bottom row and Fig. 14, right-hand column), the discrepancies from the observations are generally (though not always) smaller than in the case A 0 = 0.06 m 2 s −1 . The improved agreement in the wave vector is to be expected, but not so the agreement in the phase speed, since no effort has been made to optimize the latter. Most importantly, though, the comparison is now valid (i.e. apples-to-apples) because growing modes were found in all cases where canopy waves were observed.

Conclusions
Applying linear stability theory to boundary layer observations made over an orchard, we have explored the relationship between canopy waves and KHI. Our results support two main conclusions: -Canopy waves may be understood as KH billows only when the effects of ambient turbulence are accounted for.
-Canopy waves occur in relatively calm conditions where ambient turbulence is weak enough to permit their growth. This state is often reached just before sunrise.
Besides the inevitable measurement errors, discrepancies between theory and observations may result from three main shortcomings of the analysis.
-Questionable assumptions underlie the linear stability analyses, e.g. perturbations from parallel flow are not truly infinitesimal and the background flow is neither steady nor horizontally homogeneous. Unsteadiness may be addressed by looking for optimal transient modes of instability (Kaminski et al. 2014), while inhomogeneity requires a nonmodal approach such as large-eddy simulation (Moeng and Sullivan 2014;Patton et al. 2016). -The present analyses ignore the vertical variation of turbulence, which generally peaks at treetop height, an extension that will be implemented in the future. -Our estimates of wavelength, period and phase velocity from REAL imagery are subjective. They were a first attempt to capture the case to case variability, ignoring in particular the spatial and temporal variability that may occur within a given case. Objective methods are now under study.
Canopy waves are a significant factor in the regulation of atmospheric properties in forest canopies, especially at night or over cold surfaces, for example, at high latitudes in winter. Those properties include temperature and trace gases such as water vapor, methane, and carbon dioxide (Fitzjarrald and Moore 1990). KHI represents a more general class of phenomena which create and maintain turbulent fluxes in diverse geophysical fluid environments including oceans (Smyth and Moum 2012;Smyth 2020), estuaries (Tedford et al. 2009;Geyer et al. 2010) and the magnetospheres of Earth and other planets (Michael et al. 2021) as well as the atmosphere. The identification of canopy waves as a form of KHI suggests that, besides their intrinsic importance, forest canopies can provide a natural laboratory where KHI occurs predictably and is readily accessible for study. Topics for future study include comparison of the breaking mechanisms suggested by the asymmetric structure of observed billows with predictions from direct simulations of KHI (e.g., Espinoza-Ruiz et al. 2023;Mashayek and Peltier 2012a, b), and comparison of length scale ratios that characterize billow geometry with values found in the ocean (Smyth and Moum 2000). regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix 1: A Linear, Normal-Mode Stability Analysis of Stratified, Parallel Shear Flows Including Ambient Turbulence and Vegetation Effects
For a small-amplitude, normal mode disturbance with horizontal wave vector k evolving on a steady, parallel, horizontal background flow u h , mode evolution is determined solely by the component of the background flow parallel to the wave vector. Therefore, for a given k, we consider the two-dimensional problem of a parallel shear flow U (z) = k · u h /κ disturbed by the velocity perturbation {u , w } and the buoyancy perturbation b . Applying the Boussinesq approximation, density is decomposed as ρ = ρ 0 +ρ(z) + ρ (x, z, t), after which the buoyancy, b = g(ρ 0 − ρ)/ρ 0 , becomes b = B(z) + b , where B = −gρ/ρ 0 and b = −gρ /ρ 0 . In this context the primed quantities are treated as small perturbations.
Viscosity and diffusivity may be molecular in origin, or they may represent the effects of turbulence on spatial scales much smaller than the wave scale κ −1 . Following (Liu et al. 2012), the viscosity and diffusivity are allowed to have separate horizontal {A h , K h } and vertical {A v , K v } parts, and to vary in the vertical.
We also parameterize the transfer of momentum and buoyancy between the vegetation and the flow perturbation following (Lee 1997). The Boussinesq equations, linearized about the background flow, are then: where x is the horizontal coordinate parallel to k. The final term in each equation is the parameterized flow-vegetation flux (Lee 1997). The coefficients C d and C h are constants while plant area density L f varies with height. For the present application, the turbulence is assumed to be isotropic and to have Prandtl number unity, so that A v = K v = A h = K h . The coefficients C d and C h are either {0.15, 0.10} or zero. The vertical structure of L f is given by (6).
Perturbations u , w and b are assumed to have the normal mode form (2). The linearized normal-mode equations are then written as a generalized differential eigenvalue problem: where is the Laplacian in normal-mode form. The operators expressing viscous and diffusive effects of turbulence are: while vegetation effects are represented by: In the inviscid, nondiffusive, zero-vegetation limit reduces to the TG equation. Specification of the problem is completed by the imposition of a rigid, impermeable boundaryŵ = 0, dŵ/dz = 0 at the surface z = 0 and a frictionless boundaryŵ = d 2ŵ /dz 2 = 0 at z = 50 m. Buoyancy is fixed (b = 0) at both boundaries. The problem is solved on a stretched grid with minimum spacing 0.35 m between 10 and 30 m height and maximum 2 m in an upper sponge layer between 30 m and 50 m. Approximating z-derivatives as sixth-order finite differences, we convert the differential eigenvalue problem (10), together with the boundary conditions, to an algebraic eigenvalue problem as described by Lian et al. (2020). The eigenvalue problem is solved using standard numerical methods.