# Modelling waving crops using large-eddy simulation: Comparison with experiments and a linear stability analysis

Sylvain Dupont, Frederick Gosselin, Charlotte Py, Emmanuel de Langre, Pascal Hémon, Yves Brunet

<figure><img src="image_1.png" /></figure>

## To cite this version:

Sylvain Dupont, Frederick Gosselin, Charlotte Py, Emmanuel de Langre, Pascal Hémon, et al.. Modelling waving crops using large-eddy simulation: Comparison with experiments and a linear stability analysis. Journal of Fluid Mechanics, 2010, 652, pp.5-44. 10.1017/S0022112010000686. hal-01020671

HAL Id: hal-01020671

https://polytechnique.hal.science/hal-01020671

Submitted on 9 Jul 2014

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modelling waving crops using large-eddy simulation: comparison with experiments and a linear stability analysis

S. DUPONT<fnref n="1" />†, F. GOSSELIN<fnref n="2" />, C. PY<fnref n="3" />, E. DE LANGRE<fnref n="2" />, P. HEMON<fnref n="2" /> AND Y. BRUNET<fnref n="1" />

(Received 1 July 2009; revised 1 February 2010; accepted 1 February 2010; first published online 15 April 2010)

In order to investigate the possibility of modelling plant motion at the landscape scale, an equation for crop plant motion, forced by an instantaneous velocity field, is introduced in a large-eddy simulation (LES) airflow model, previously validated over homogeneous and heterogeneous canopies. The canopy is simply represented as a poroelastic continuous medium, which is similar in its discrete form to an infinite row of identical oscillating stems. Only one linear mode of plant vibration is considered. Two-way coupling between plant motion and the wind flow is insured through the drag force term. The coupled model is validated on the basis of a comparison with measured movements of an alfalfa crop canopy. It is also compared with the outputs of a linear stability analysis. The model is shown to reproduce the well-known phenomenon of ‘honami’ which is typical of wave-like crop motions on windy days. The wavelength of the main coherent waving patches, extracted using a bi-orthogonal decomposition (BOD) of the crop velocity fields, is in agreement with that deduced from video recordings. The main spatial and temporal characteristics of these waving patches exhibit the same variation with mean wind velocity as that observed with the measurements. However they differ from the coherent eddy structures of the wind flow at canopy top, so that coherent waving patches cannot be seen as direct signatures of coherent eddy structures. Finally, it is shown that the impact of crop motion on the wind dynamics is negligible for current wind speed values. No lock-in mechanism of coherent eddy structures on plant motion is observed, in contradiction with the linear stability analysis. This discrepancy may be attributed to the presence of a nonlinear saturation mechanism in LES.

## 1. Introduction

Turbulent wind flows over vegetation canopies are dominated by intermittent, energetic downward-moving gusts, known as large coherent eddy structures, which evolve within unorganized random background turbulence. It has been demonstrated that these coherent structures are generated by processes similar to those occurring in a plane-mixing layer flow (Raupach, Finnigan & Brunet 1996). In response to the forcing of these structures, plants sway like damped harmonic oscillators

<footnote n="1">¹INRA, UR1263 EPHYSE, 71 avenue Edouard Bourlaux, F-33140 Villenave d'Ornon, France</footnote>

<footnote n="2">²Département de Mécanique, LadHyX, CNRS-Ecole Polytechnique, F-91128 Palaiseau, France</footnote>

<footnote n="3">³MSC, UMR 7057 CNRS-Université Paris-Diderot, F-75205 Paris cedex 13, France</footnote>

<footnote>† Email address for correspondence: sdupont@bordeaux.inra.fr</footnote>
<figure><img src="image_3.png" /><figcaption>FIGURE 1. (a) Visualization of coherent waving patches (white patches) on an alfalfa field under wind forcing. (b) Experimental set-up for the measurement of the wind-induced motion of an alfalfa crop performed by Py et al. (2005). (c) Velocity field of the upper surface of the crop deduced from video-recording. Figure adapted from Py et al. (2005).</figcaption></figure>

excited by intermittent impulsive loading (Finnigan 1979; Gardiner 1995). A striking visualization of the interaction between coherent structures and plant movements on windy days is known as ‘honami’, which refers to wave-like crop motions (Inoue 1955; Finnigan & Mulhearn 1978a). One example of ‘honami’ waves was recorded by Py et al. (2005) over an alfalfa crop canopy (figure 1a). While turbulent wind flow over plant canopies has been widely studied (see Finnigan 2000 for a review), coherent motions of the canopies themselves and their strong interactions with the wind flow have received little attention so far.
The main motivations for better understanding wind-induced plant motion in a fully-coupled way were reported in the recent review of de Langre (2008) as (i) environmental applications such as management practices aimed at limiting forest and crop damage due to windthrow and lodging in strong wind conditions, (ii) agronomic applications aimed at improving plant growth and biomass production through the thigmomorphogenesis effect (Jaffe 1973; Moulia & Combes 2004) and (iii) image synthesis applications. Wind-induced plant movements are particularly complex since agricultural landscapes often exhibit large spatial variability caused by the presence of adjacent crop fields, clearings, roads, forest patches of various height, etc. Such heterogeneities exert influence on the turbulent fields in the lower atmosphere, and thereby on canopy motion. Because of the complexity of the various processes responsible for plant motion in heterogeneous landscapes, modelling both plant and flow dynamics appears necessary for quantifying plant motion according to the position in the landscape. Developing and validating such a computational tool for better understanding wind-plant interaction is the subject of the present paper.

## 1.1. *Wind flow and plant motion models*

Airflow within and above vegetation canopies has been investigated through several Reynolds-averaged type models (e.g. Li & Lin 1990; Green 1992; Liu et al. 1996; Foudhil, Brunet & Caltagirone 2005; Dupont & Brunet 2006). However these models only simulate mean fields, which is a serious limitation for studying wind-induced plant motion, that is generated primarily by wind gusts. Large eddy simulation (LES) techniques allow one to analyse canopy turbulence in much greater details since eddy motions larger than twice the grid mesh are explicitly solved and only subgrid-scale (SGS) eddy motions are modelled. Provided that the grid is fine enough, LES therefore allows one to have access to instantaneous dynamic fields and consequently is capable of reproducing wind gusts in a plant canopy. Over the last decade it has been demonstrated that the LES technique reproduces the main features of turbulent flow observed over homogeneous vegetation canopies (Shaw & Schumann 1992; Kanda & Hino 1994; Dwyer, Patton & Shaw 1997; Shen & Leclerc 1997; Su et al. 1998; Su, Shaw & Paw, U 2000; Watanabe 2004; Dupont & Brunet 2008c), downwind from forest leading edges (Yang et al. 2006a,b; Dupont & Brunet 2008a,b, 2009) as well as over forested hills (Dupont, Brunet & Finnigan 2008; Ross 2008). However, in these airflow models, the canopy is usually represented by a simple drag force term in the momentum equation, without accounting for plant motion; in other words, plant elements are considered smaller than the airflow grid cells and fixed in space. In the canopy-flow literature it is indeed usually considered that plant motion has negligible impact on the flow both within and above the canopy. The main advocated reason is that turbulent structures induced by plant motion are much smaller than the main coherent structures of the canopy flow, as the former scale with the dimension of plant elements while the latter scale with canopy height.

A large range of models have also been developed to simulate plant motion, from simple mass-spring-damper models (Finnigan & Mulhearn 1978b; Mayer 1987; Flesch & Grant 1992; Flesch & Wilson 1999; Farquhar, Wood & van Beem 2000; Doaré, Moulia & De Langre 2004; Py, de Langre & Moulia 2004, 2006; Gosselin & de Langre 2009) to complex dynamic models based on the finite element method (Kerzenmacher & Gardiner 1998; Ikeda, Yamada & Toda 2001; Sellier, Fourcaud & Lac 2006; Sellier, Brunet & Fourcaud 2008; Rodriguez, de Langre & Moulia 2008). In the simplest models plants are represented as one- or two-dimensional oscillating rods, and for the most complex ones by flexible beams with branches. Complex models are
usually applied at the individual plant scale, for simulating tree motion for example, whereas simpler models are applied at canopy scale for simulating crop motion. In the latter case the canopies are considered as poroelastic systems since the plants are not individualized (de Langre 2008). Plant or crop models may be forced by analytical functions (wind pulse, sinusoidal wavetrain), or by time series of measured wind velocity. Feedback from plant motion to the wind flow is not usually considered except in a few models like those of Finnigan & Mulhearn (1978b), Ikeda et al. (2001), Py et al. (2006) and Gosselin & de Langre (2009), which are discussed in the next subsection. This review would not be complete without mentioning mechanistic models such as HWIND (Peltola et al. 1999), GALES (Gardiner, Peltola & Kellomaki 2000) or FOREOLE (Ancelin, Courbaud & Fourcaud 2004), which have been developed to quantify tree stability to windload, in forest management perspectives. Unlike the previous ones, these models are static and turbulence is only accounted for through a gust factor deduced from wind-tunnel measurements.

## 1.2. Interaction between wind flow and plant motion

As already mentioned, a few studies only have focused on the interactions between wind flow and plant motion. Firstly, Finnigan & Mulhearn (1978a) obtained qualitative results using a flexible canopy model of a wheat field, placed in a wind tunnel. They observed a peak in the velocity spectra at the waving frequency of individual plant models, which imply that plant motions may modify the wind flow. Finnigan & Mulhearn (1978b) then developed a mathematical model based on a linearized one-dimensional momentum equation and a plant motion equation coupled through a drag force. With this model they investigated the response of the wind flow fluctuations to varying excitation frequency in order to study the impact of plant motion on wind-flow fluctuations within the canopy, and how it varies with plant spacing, flexibility, leafiness and mean wind speed. Their study confirmed that plant motion in a dense flexible canopy model may alter the wind flow. More recently, Ikeda et al. (2001) introduced in a LES model a motion equation for a canopy of flexible plants, in order to analyse the impact of plant motion on turbulence over a reed field. They observed that the periodicity of wind vortex generation is reduced when the plants are assumed flexible. However their study only presents qualitative information and use a questionable ‘two-dimensional LES’ model. More recently, Ghisalberti & Nepf (2006) studied the structure of coherent eddies over rigid and flexible canopies in a flume tunnel.

Possibly the most detailed dataset available to date on crop motion comes from the video-recording experiment performed by Py et al. (2006) over alfalfa and wheat fields, which allowed them to characterize the spatio-temporal movements of crops subjected to wind. Their experimental technique, based on image correlation analysis of crop motion, is described in Py et al. (2005). Py et al. (2006) completed this experiment with a linear stability analysis performed with a two-dimensional analytical model, fully coupling a mixing-layer flow with an oscillating vegetation canopy through a drag force. Canopy motion was driven by the Kelvin-Helmholtz instability of the modelled flow, instead of being forced by an imposed oscillating flow as in Finnigan & Mulhearn (1978b). With their linear stability analysis, Py et al. (2006) observed a lock-in mechanism similar in form to what is observed in vortex-induced vibrations. As the mean wind speed increases, the frequency of the Kelvin-Helmholtz instability increases; it deviates as it approaches the plant frequency and locks onto it. Hence, within a specific range of wind velocity, the flow and the vegetation canopy move in phase. This lock-in mechanism was further studied by Gosselin & de Langre (2009) on aquatic plants in a water stream, using a revisited version of the model
of Py et al. (2006) and the data of Ghisalberti & Nepf (2006). Py et al. (2006) showed that with a reasonable mean wind profile, their stability analysis over alfalfa and wheat crops can predict the wavelengths deduced from the video recordings of Py et al. (2005). However, the broken-line mean wind velocity profile imposed in Py et al. (2006) is a rough approximation of reality, with arbitrary values attributed to its constitutive parameters. This limitation, combined with the fact that no measurements led to experimental points outside the lock-in range, leaves the existence of lock-in unconfirmed. It therefore follows from these various studies that the possible impact of plant motion on turbulence remains unclear.

### 1.3. Objectives

The first goal of the present study is to introduce a novel three-dimensional computational model strongly coupling plant motion and turbulent wind flow. For this purpose, an equation for plant motion was introduced into an atmospheric LES model. The canopy is represented as a poroelastic continuum. This representation is similar in its discrete form to an infinite row of identical mechanical oscillators where only a linear mode of plant vibration is considered. The two-way coupling between plant motion and wind flow occurs through the drag force term. The model is validated against video recordings of alfalfa crop motion performed by Py et al. (2006).

The second goal of this study is to understand, using the LES model, the mechanisms governing ‘honami’ and whether plant motion influence wind turbulence. In particular, we would like to address the two following questions. Firstly, can coherent crop motions be considered as direct signatures of coherent eddies? In other words, is it possible to deduce spatial and temporal information on coherent eddies, or simply information on the wind flow, from crop motion modelling or video recordings? Secondly, can crop motion alter wind flow? More specifically, can a lock-in mechanism occur in LES simulations? This leads us to investigate the differences between three-dimensional, nonlinear effects simulated by LES from the linear two-dimensional interaction mechanisms accounted for by the stability analysis.

The modified LES model is first presented in § 2, along with the numerical simulation set-up, the waving crop experiment of Py et al. (2005, 2006) used to validate the LES model, and the equations of the linear stability analysis used for comparison with the LES model. We then analyse, and test against measurements, the magnitude and velocity of plant deflection over a range of wind velocity at canopy top (§ 3), as well as the spatial and temporal characteristics of the main coherent crop motions (§ 4). We further investigate the interaction between organized crop motion and coherent eddy structures by comparing their main properties (§ 5) and by looking at the potential impact of crop motion on the turbulent wind flow (§ 6). The lock-in mechanism is discussed in § 7 by comparing LES with a linear stability analysis, and we conclude in § 8.

## 2. Method

### 2.1. *The waving crop experiment*

The measurements performed by Py et al. (2005) in a field of waving alfalfa (Medicago sativa L. cv Mercedes) are used in the present study to validate our model. We only present here the main characteristics of their experiment. Full details can be found in Py et al. (2005).

Canopy motion was video-recorded at a frequency of 25 frames per second in various wind conditions (figure 1b), during sequences of 10–30 s. Throughout each recorded sequence wind velocity was measured with a hot-wire anemometer located just above the crop surface. After correction of the images from distortion,
<table><caption>TABLE 1. Structural and mechanical properties of the plants considered in this study.</caption><thead><tr><th>Parameter</th><th>Designation</th><th>Alfalfa</th><th>Wheat</th></tr></thead><tbody><tr><td>C<sup>canopy</sup><sub>d</sub></td><td>Canopy drag coefficient</td><td>0.20</td><td>0.20</td></tr><tr><td>f<sub>0</sub></td><td>Natural vibration frequency (Hz)</td><td>1.05</td><td>2.50</td></tr><tr><td>h</td><td>Mean canopy height (m)</td><td>0.69</td><td>0.69</td></tr><tr><td>l</td><td>Mean plant spacing (m)</td><td>0.05</td><td>0.05</td></tr><tr><td>LAI</td><td>Leaf-area index</td><td>3.00</td><td>3.00</td></tr><tr><td>m</td><td>Plant mass (kg)</td><td>0.014</td><td>0.007</td></tr><tr><td>ξ</td><td>Damping coefficient</td><td>0.0875</td><td>0.0859</td></tr></tbody></table>

a two-dimensional spatio-temporal horizontal velocity field of the crop surface was deduced by a correlation analysis based on standard particle image velocimetry (PIV) algorithms, where small-scale plant heterogeneities play the role of natural tracers (figure 1c). The main characteristics of the coherent structures observed in canopy motion (spatial wavelength and temporal frequency) were deduced from bi-orthogonal decompositions (BOD) of the crop velocity field at various wind speeds. Additionally, the structural and mechanical properties of six individual plants were measured. Their mean values are reported in table 1.

The atmosphere stability was unfortunately not measured during this experiment. However, since measurements were performed during daytime, the atmosphere should not be stable; and since the magnitudes of the wind velocity at canopy top were mostly larger than 1 m s⁻¹, free convective conditions should not be present. Therefore, we think that it is reasonable, at the scale of the crop canopy and for studying plant motion, to consider hereafter a neutral stratification of the atmosphere for the validation of our model against the present measurements.

## 2.2. Large-eddy simulation model equations

In order to simulate the wind flow and canopy plant motion, we use the Advanced Regional Prediction System (ARPS, version 5.1.5) originally developed at the Center for Analysis and Prediction of Storms (CAPS), University of Oklahoma, for the explicit prediction of convective and cold-season storms as well as weather systems. A detailed description of the standard version of the model and its validation cases are available in the ARPS User's Manual (Xue et al. 1995) and in Xue, Droegemeier & Wong (2000) and Xue et al. (2001).

ARPS is a three-dimensional nonhydrostatic compressible model where Navier-Stokes equations are written in the so-called Gal-Chen, or terrain-following coordinates. The grid is orthogonal in the horizontal direction and stretched in the vertical. The model solves the conservation equations for the three wind velocity components, pressure, potential temperature and water (water vapour, cloud water, rainwater, cloud ice, snow and graupel). Wind components and atmospheric state variables (air density, pressure and potential temperature) are split into a base state (hereafter represented by over-barred variables) and a deviation (double-primed variables). The base state is assumed horizontally homogeneous, time invariant and hydrostatically balanced. At high spatial resolution the conservation equations are implicitly filtered towards the grid, in order to separate the small scales from the large scales.

The large-eddy simulation (LES) approach used by ARPS consists in resolving explicitly all turbulent structures larger than the filter scale, while smaller turbulent structures, i.e. SGS turbulent motions, are modelled using an eddy viscosity approach.
As this type of model is entirely dissipative, it does not account for energy backscattering from small to large scales. The eddy viscosity is represented as the product of a length scale and a velocity scale characterizing the SGS turbulent eddies. An SGS turbulent kinetic energy (TKE) conservation equation is solved so as to obtain the representative velocity scale, and the length scale is derived from the grid spacing. Different horizontal and vertical grid spacings lead to using different horizontal and vertical length scales. The SGS model is described in Appendix A.

Recently, Dupont & Brunet (2008c) modified the model so as to simulate turbulent flows at very fine scales (0.1h, where h is the mean canopy height) within canopies made of fixed plants. The mean turbulent fields and the development of coherent structures, as simulated by this modified version of ARPS, were successfully validated against field and wind-tunnel measurements, over homogeneous canopies (Dupont & Brunet 2008c), a simple forest-clearing-forest pattern (Dupont & Brunet 2008a,b, 2009) and a forested hill (Dupont et al. 2008). Here, the model is further extended so as to simulate plant motion and its interaction with the wind flow.

Within the vegetated layer, the shear stress at canopy top generates eddies larger than the eddies formed in the wake of the vegetation elements, and TKE dissipates through the smallest eddies (Kolmogorov scale). Shear-type structures are explicitly resolved by the model while wake-type structures are modelled. To account for the presence of vegetation on the wind flow, a drag-force approach was implemented by adding a pressure and viscous drag force term in the momentum equation (2.1), and by adding a sink term in the equation for SGS TKE (see (A 6)) in order to represent the acceleration of the dissipation of turbulent eddies in the inertial subrange. As all simulations in this study were performed in a dry neutrally stratified flow over a flat terrain, the momentum equation presented hereafter is written in Cartesian coordinates for a dry atmosphere. Although the atmosphere is assumed neutral, the potential temperature equation (not shown) has to be solved because turbulent motions are activated through initial turbulent perturbations. The momentum equation, written for a Boussinesq fluid and using the Einstein summation convention, therefore reads

$$ \frac{\partial \tilde{u}_i}{\partial t} + \tilde{u}_j \frac{\partial \tilde{u}_i}{\partial x_j} = -\frac{1}{\bar{\rho}} \frac{\partial}{\partial x_i} \left( \tilde{p}'' - \alpha_{div} \frac{\partial \bar{\rho} \tilde{u}_j}{\partial x_j} \right) - g \left( \frac{\tilde{\theta}''}{\bar{\theta}} - \frac{c_p}{c_v} \frac{\tilde{p}''}{\bar{p}} \right) \delta_{i3} \\ - \frac{\partial \tau_{ij}}{\partial x_j} - \frac{C_D}{l^2} \left| \tilde{u}_i - (1 - \delta_{i3}) \frac{z}{h} \frac{\partial \tilde{q}_i}{\partial t} \right| \left( \tilde{u}_i - (1 - \delta_{i3}) \frac{z}{h} \frac{\partial \tilde{q}_i}{\partial t} \right), \quad (2.1) $$

where the overtilde indicates the filtered variables or grid volume-averaged variables. Here, t is time; $x_i$ ($x_1=x$, $x_2=y$, $x_3=z$) refer to the streamwise, lateral and vertical directions, respectively; $u_i$ ($u_1=u$, $u_2=v$, $u_3=w$) is instantaneous velocity component along $x_i$; $q_i$ ($q_1=q_x$, $q_2=q_y$) is instantaneous plant displacement component along $x_i$ at canopy top; $\delta_{ij}$ is Kronecker delta, $\alpha_{div}$ is a damping coefficient meant to attenuate acoustic waves; p is air pressure; $\rho$ is air density; g is acceleration due to gravity; $\theta$ is potential temperature; $\tau_{ij}$ is subgrid stress tensor defined in Appendix A; $c_p$ and $c_v$ are specific heat of air at constant pressure and volume, respectively.

The terms on the right-hand side of (2.1) represent, respectively, the pressure-gradient force term, the buoyancy term, the turbulent transport term, and the drag force term induced by the vegetation. The latter term is proportional to the relative velocity between the wind $\tilde{u}_i$ and the plant deflection velocity ($z/h \partial \tilde{q}_i / \partial t$). In this term, l is the average plant spacing, $C_D = C_d^{canopy} A_f^{plant}$, where $C_d^{canopy}$ and $A_f^{plant}$ are
<figure><img src="image_4.png" /><figcaption>FIGURE 2. (a) Schematic representation of crop plants as oscillating stems under wind forcing where $q_i$ is the plant displacement at canopy top in direction $i$. (b) Vertical profile of the frontal area density of the alfalfa canopy considered in the simulations.</figcaption></figure>

the mean canopy drag coefficient and the mean plant frontal area density ($m^2 m^{-1}$), respectively.

Plants in a crop canopy can be seen as identical mechanical oscillating stems with two degrees of freedom. This was shown by Finnigan & Mulhearn (1978b) and confirmed by vibration tests performed by Py et al. (2006) on alfalfa and wheat crops. Following the modal analysis, the deformation of plant stems can be decomposed into a set of various vibration modes so that its displacement $q_i$ in direction $i$ is the sum of the contributions of each vibration mode: $q_i(t) = \sum_n \lambda_i^n(t) \varphi_i^n$, where $\varphi_i^n$ represents the mode shape $n$ of the stem and $\lambda_i^n$ its associated displacement. As crop plants, such as wheat and alfalfa, have a very slender shape, most of their dynamics may be represented by the first mode of vibration. Hence, only the fundamental mode of plant vibration is considered in this study, and it is further approximated by a linear mode shape, $\varphi_i = z/h$, consistently with the representation used by Doaré et al. (2004) and Py et al. (2006) for alfalfa and wheat plants and by Flesch & Grant (1991) for corn plants. The use of a linear mode shape means that the plant deformation is represented by a function varying linearly with $z$. Hence, the angular displacement of the stem is constant along its height (see figure 2a). Doaré et al. (2004) observed that this linear mode shape is well suited for modelling alfalfa plant motion. Note that the shape of the mode only plays a role in weighting in space the coupling between the oscillating canopy and the flow. With this approach, plant stems are only characterized by their height $h$, mass $m$, non-dimensional damping coefficient $\xi$ and natural vibration frequency $f_0$. Furthermore, the angular displacements of plants $\phi_i$ in direction $i$ are assumed sufficiently small to consider that $\sin \phi_i \approx \phi_i$ and $\cos \phi_i \approx 1$. Hence, the kinematics of an isolated flexible stem under wind load is described by the following simple mass-spring-damper equation:

$$ \underbrace{\left(\frac{1}{3}mh^2\right)}_{(i)} \frac{\partial^2 \phi_i}{\partial t^2} + \underbrace{c \frac{\partial \phi_i}{\partial t}}_{(ii)} + \underbrace{r \phi_i}_{(iii)} - \underbrace{\left(\frac{1}{2}mgh\right)}_{(iv)} \phi_i = \underbrace{\bar{\rho} \int_0^h C_D \left|\tilde{u}_i - z \frac{\partial \phi_i}{\partial t}\right| \left(\tilde{u}_i - z \frac{\partial \phi_i}{\partial t}\right) h \phi_i dz}_{(v)} \quad (2.2) $$
where $i \in \{1, 2\}$, c and r are plant damping and stiffness coefficients, respectively. The terms on the left-hand side of (2.2) represent, respectively, the inertia (i), damping (ii), stiffness (iii) and gravity (iv) terms. The term on the right-hand side is the wind-induced drag term (v); here, the mean plant drag coefficient is assumed equal to that of the canopy.

At crop scale, the canopy can be seen as a succession of infinite rows of identical bi-dimensional mechanical oscillating rigid stems, where each stem displacement is solved from (2.2). In order to use the same horizontal resolution between wind flow and plant motion, the canopy is not seen as a succession of individual stems but as a poroelastic continuum medium whose motion is described by the grid volume-averaged displacement $\tilde{q}_i(x, y, t)$ at canopy top. With some notation simplifications, the continuous form of (2.2) writes as follows:

$$M \frac{\partial^2 \tilde{q}_i}{\partial t^2} + C \frac{\partial \tilde{q}_i}{\partial t} + R \tilde{q}_i = \bar{\rho} \int_0^h C_D \left| \tilde{u}_i - \frac{z}{h} \frac{\partial \tilde{q}_i}{\partial t} \right| \left( \tilde{u}_i - \frac{z}{h} \frac{\partial \tilde{q}_i}{\partial t} \right) \varphi_i \, dz, \quad (2.3)$$

where $i \in \{1, 2\}$, $q_i = h\phi_i$, $M = m/3$, $C = c/h^2$ and $R = r/h^2 - mg/(2h)$. Here, M is the mass and C and R are, respectively, the damping and stiffness coefficients of plant. The damping coefficient is computed from $c = 4\pi m h^2 f_0 \xi / 3$ and the stiffness coefficient is deduced from the relationship $f_0 = R/(4\pi^2 M)$, which leads to $r = 4\pi^2 m h^2 f_0^2 / 3 + mg/2$. Compared to the oscillator equation used by Py et al. (2006) and presented in §2.4 (see (2.7)), the present equation is two-dimensional and the drag term is not linearized.

Interactions between neighbouring plants are neglected as they were found negligible by Py et al. (2006) for an alfalfa canopy. Nevertheless, with this continuous form of the equation for crop motion, elastic contacts between plants could be easily considered in the future, through an additional term depending on the second spatial derivative of plant displacement, transforming equation (2.3) into a wave-like equation (Doaré et al. 2004).

As already mentioned, for this initial version of the model, only small plant displacements are considered. This assumption allows to consider that the motion of a plant always occurs inside the grid box of its rest position. Hence, the wind velocity $\tilde{u}_i$ appearing in the drag term of (2.3) is simply the grid volume-averaged velocity solved from (2.1) where is located the plant element at rest. No interpolation of the wind velocity within the grid cell to the position of the plant is considered as well as no account for SGS velocity. This assumption appears reasonable for horizontal grid size larger than horizontal plant displacements.

As reviewed by de Langre (2008), the deformation shape of plant through wind load becomes more streamlined, reducing the frontal area density and the pressure component of drag, and so affecting the drag load. This effect is usually accounted by modifying the dependence of the drag load with the velocity from a square to a linear dependence with increasing wind speed. In our model, as only small plant displacements are considered, the square dependence of the drag term is always considered and the vertical profile of $C_D$ is taken constant with wind velocity and plant deflection.

To summarize, our crop-plant motion model considers the following assumptions: (i) identical properties for all plants of a crop, (ii) linear mode shape of plant deformation, (iii) small plant deflections, (iv) plant motion inside the same fluid grid cell during the simulation, (v) no interaction between neighbouring plants and (vi) no streamlining effect due to plant deformation with increasing wind velocity. The
last three assumptions are a consequence of the third assumption, i.e. small plant
deflections.

For the sake of clarity the overtilde on $\tilde{u}_i$ and $\tilde{q}_i$ will be omitted from now on, and
the plant displacement velocity will be noted as $\zeta_i = \partial q_i / \partial t$.

2.3. Numerical details of large-eddy simulations

Four three-dimensional simulations were performed over a homogeneous continuous
alfalfa crop canopy with different values of canopy-top wind speed $U_h$ ranging from 1
to 4 m s⁻¹. Such velocity values are currently observed over crop canopies, and extreme
values of wind speed were not considered in this study. These four simulations are
hereafter referred to as Cases 1–4. Properties of the alfalfa plant were chosen as
similar to those of Py *et al.* (2006) (table 1). Plant height $h$ was set to 0.69 m and
the average plant spacing $l$ to 0.05 m. The vertical distribution of the frontal area
density $A_f^{\text{plant}}$ was assumed identical to the average vertical mass distribution of the
six plants measured by Py *et al.* (2006). This leads to a constant profile of $A_f^{\text{plant}}$
within the lower canopy and a linear decrease in the upper part (see figure 2b). The
magnitude of $A_f^{\text{plant}}$ was chosen so as to provide a leaf-area index (LAI) of about 3
($\int_0^h A_f^{\text{plant}} l^{-2} dz = 3$), which is typical of alfalfa crops (Russell, Marshall & Jarvis 1990).
The drag coefficient $C_d^{\text{canopy}}$ was taken as 0.2.

All simulations were performed within a unique computational domain, extending
over 30 × 15 × 8 m³. This corresponds to 200 × 100 × 65 grid points in the x, y and
z directions, respectively, and to a horizontal resolution of 0.15 m. The vertical grid
resolution is 0.08 m below 3.5 m, and above the grid is stretched using a hyperbolic
tangent function with a vertical resolution of 0.4 m at the top of the domain. This
choice of size and resolution of our computational domain results from a compromise
between constraints related to the available computational time and the spatial and
lifetime resolutions of the main eddies responsible for plant motions. These latter
turbulent structures are induced by the canopy-top mean wind shear, their horizontal
and vertical extents are of the order of $h$ and $h/3$, respectively (Finnigan 2000).
Consequently, the resolution of our domain should be sufficient for simulating such
structures. On the other hand, the size limitation of the domain does not allow
mesoscale structures or large atmospheric surface layer to be resolved since they
have a much larger spatial scales than our domain and much larger lifetime than the
duration of our simulation. The small size of our domain compared to the planetary
boundary layer should not be a problem since the atmosphere is taken as neutral
and since canopy turbulent structures are the main structures of interest in plant
motions. Hence, in our simulations, mesoscale structures should be considered as a
background average flow dynamics.

The choice of the present resolution induces a grid aspect ratio ($\Delta x/\Delta z$, $\Delta y/\Delta z$) of 1.9 at the bottom of the domain and 2.7 at its top, which leads to use the 1.5-order closure scheme in its anisotropic form where two mixing lengths are computed for horizontal and vertical turbulent diffusion (see Appendix A). As reported by Chow *et al.* (2006), a too large aspect ratio may induce numerical errors in the horizontal gradient terms (Mahrer 1984) as well as some distortion of the resolved eddies, in particular, over complex terrain (Kravchenko, Moin & Moser 1996). Our aspect ratio of 1.9 close to the surface appears much smaller than the value of 10 used with ARPS by Chow *et al.* (2006), Weigel *et al.* (2006) and Weigel, Chow & Rotach (2007) over a complex terrain and the value of 4.7 used by Cassiani, Katul & Albertson (2008) over an heterogeneous canopy in a flat terrain. Furthermore, these authors found no significant differences with simulations performed with a smaller aspect ratio. From
this, we believe that in the present study our aspect ratio is quite reasonable and sufficient to obtain a realistic estimate of turbulence structures over the crop canopy.

The lateral boundary conditions are periodic, the bottom boundaries are treated as rigid and the surface momentum flux is parameterized by using bulk aerodynamic drag laws. A 2.5 m deep Rayleigh damping layer is used at the upper boundary in order to absorb upward-propagating wave disturbances and to eliminate wave reflection at the top of the domain. Additionally, the flow is driven by a depth-constant geostrophic wind corresponding to a base-state wind at the upper boundary. The velocity fields are initialized using a meteorological pre-processor (Pénelon, Calmet & Mironov 2001) with a constant vertical profile of potential temperature and a dry atmosphere. The plant motion equation was resolved from an explicit time integration scheme using the same time step of 0.0015 s as the momentum equation. This time step, which is much smaller than the period of natural vibration of alfalfa plants (0.95 s), should insure a natural plant swaying.

After the flow has reached an equilibrium state, wind and plant turbulent statistics were computed from a horizontal- and time-averaging procedure. Horizontal averaging was performed over all x and y locations at each considered z, and time averaging was performed over 300 samples collected during a 30 s period. Consequently, wind velocity components as well as plant displacement and velocity components can be decomposed into $\phi_i = \langle \phi_i \rangle_{xyt} + \phi'$, where $\phi_i$ is either $u_i$, $q_i$ or $\zeta_i$, the symbol $\langle \rangle_{xyt}$ denotes the time and space average and the prime denotes the deviation from the averaged value.

## 2.4. Linear stability analysis equations

In order to emphasize the basic mechanisms that govern the complex flow modelled in the nonlinear three-dimensional LES, we also performed a linear stability analysis with the two-dimensional model of Py et al. (2006) which couples a mixing-layer flow with an oscillating canopy. For the sake of realism and to perform easier comparisons, the model of Py et al. (2006) was modified so as to include the effects of eddy-viscosity and use a more realistic continuous mean wind profile than a broken-line velocity profile.

We thus study the linear stability of the base flow (overbarred variables) to small perturbations $u'(x, z, t)$, $w'(x, z, t)$, $p'(x, z, t)$ and $q'_x(x, t)$, respectively, the x and z components of the perturbation velocity, the perturbation pressure and the plant displacement at canopy top. The base flow is characterized by the mean wind profile $\bar{U}(z)$ and the isotropic eddy viscosity $\bar{\nu}_t(z)$, which only depend on z and are imposed on the system. Perturbations in $\nu_t$ are neglected for simplicity. The linearized equations for momentum conservation in x and z of Py et al. (2006), with the added terms of eddy viscosity, read as follows, along with the continuity equation and the linearized oscillator equation governing the dynamics of the canopy:

$$ \frac{\partial u'}{\partial t} + \bar{U} \frac{\partial u'}{\partial x} + w' \frac{\partial \bar{U}}{\partial z} = -\frac{\partial p'}{\partial x} + \bar{\nu}_t \left( \frac{\partial^2 u'}{\partial x^2} + \frac{\partial^2 u'}{\partial z^2} \right) \\ + \frac{\partial \bar{\nu}_t}{\partial z} \left( \frac{\partial u'}{\partial z} + \frac{\partial w'}{\partial x} \right) - 2 \frac{C_D}{l^2} \bar{U} \left( u' - \frac{z}{h} \frac{\partial q'_x}{\partial t} \right), \quad (2.4) $$
$$ \frac{\partial w'}{\partial t} + \bar{U} \frac{\partial w'}{\partial x} = -\frac{\partial p'}{\partial z} + \bar{v}_t \left( \frac{\partial^2 w'}{\partial x^2} + \frac{\partial^2 w'}{\partial z^2} \right) + 2 \frac{\partial \bar{v}_t}{\partial z} \frac{\partial w'}{\partial z}, \quad (2.5) $$

$$ \frac{\partial u'}{\partial x} + \frac{\partial w'}{\partial z} = 0, \quad (2.6) $$

$$ M \frac{\partial^2 q_x'}{\partial t^2} + C \frac{\partial q_x'}{\partial t} + R q_x' = 2\bar{\rho} \int_0^h C_D \bar{U} \left( u' - \frac{z}{h} \frac{\partial q_x'}{\partial t} \right) \frac{z}{h} dz. \quad (2.7) $$

Boundary conditions of no penetration and free slip at the ground ($z=0$) and at the top of the domain ($z=H$) are applied to the flow field, i.e.

$$ \left. \begin{array}{l} w' \bigg|_{z=0} = 0, \quad \left[ \frac{\partial u'}{\partial z} + \frac{\partial w'}{\partial x} \right]_{z=0} = 0, \\ w' \bigg|_{z=H} = 0, \quad \left[ \frac{\partial u'}{\partial z} + \frac{\partial w'}{\partial x} \right]_{z=H} = 0. \end{array} \right\} \quad (2.8) $$

We seek a solution to (2.4)–(2.8) in the form of a travelling wave

$$ \langle u', w', q_x' \rangle = \langle \hat{u}, \hat{w}, \hat{q}_x \rangle e^{i(kx - \omega t)} + \text{c.c.}, \quad (2.9) $$

where $k$ and $\omega$ are the streamwise wavenumber and the complex frequency and where c.c. stands for complex conjugate. Upon substitution of the travelling wave solution, (2.4)–(2.6) can be combined into a single differential equation of $\hat{w}$ and $\hat{q}_x$:

$$ \omega \left( ik^2 \hat{w} - i \frac{d^2 \hat{w}}{dz^2} \right) + ik \bar{U} \left( -k^2 \hat{w} + \frac{d^2 \hat{w}}{dz^2} \right) - ik \frac{d^2 \bar{U}}{dz^2} \hat{w} \\ - \bar{v}_t \left( k^4 \hat{w} - 2k^2 \frac{d^2 \hat{w}}{dz^2} + \frac{d^4 \hat{w}}{dz^4} \right) - 2 \frac{d \bar{v}_t}{dz} \left( \frac{d^3 \hat{w}}{dz^3} - k^2 \frac{d \hat{w}}{dz} \right) - \frac{d^2 \bar{v}_t}{dz^2} \left( \frac{d^2 \hat{w}}{dz^2} + k^2 \hat{w} \right) \\ + \frac{2}{l^2} \frac{d \hat{w}}{dz} \frac{d}{dz} (C_D \bar{U}) + \frac{2}{l^2} (C_D \bar{U}) \frac{d^2 \hat{w}}{dz^2} + \hat{q}_x \frac{2\omega k}{hl^2} \frac{d}{dz} (C_D \bar{U} z) = 0. \quad (2.10) $$

We introduce the new quantity $\zeta_x' = \partial q_x' / \partial t$ with $\zeta_x' = \hat{\zeta}_x e^{i(kx - \omega t)} + \text{c.c.}$, such that the second-order differential equation (2.7) can be written as two first-order equations:

$$ -i\omega M \hat{\zeta}_x + C \hat{\zeta}_x + R \hat{q}_x = 2\bar{\rho} \int_0^h C_D \bar{U} \left( \frac{i}{k} \frac{d\hat{w}}{dz} + i\omega \frac{z}{h} \hat{q}_x \right) \frac{z}{h} dz, \quad (2.11) $$

$$ \hat{\zeta}_x = -i\omega \hat{q}_x. \quad (2.12) $$

The $z$-function of the vertical perturbation flow velocity $\hat{w}(z)$ is discretized over the domain $[0, H]$ at $N+2$ nodes. Upon substituting a centred finite-difference scheme, the system of equations (2.12), (2.11) and (2.10) can be formulated as an eigenvalue problem:

$$ (A - \omega B) w = 0, \quad (2.13) $$

where $w = \langle \hat{\zeta}_x, \hat{q}_x, \hat{w}_1, \hat{w}_2, \dots, \hat{w}_N \rangle^\mathrm{T}$ and where the linear operators $A$ and $B$ are given in Appendix B. For a given value of $k$, (2.13) is solved for its eigenvalues. Each combination of $k$ and $\omega$ satisfying the governing equations corresponds to a mode of the system. For each mode, the complex frequency has a real and an imaginary part, $\omega = \omega_r + i\omega_i$. The real part is the oscillation frequency and the imaginary part is the temporal growth rate. If $\omega_i > 0$, the vibration mode is unstable and a small perturbation will increase exponentially. On the other hand, if $\omega_i < 0$ the vibration mode is stable and a small perturbation will decay. If $\omega_i = 0$ the mode is neutrally stable.
Once the eigenfrequencies of the system are found, we use the corresponding eigenfunctions to compute the distribution of perturbation energy. Similarly to Gosselin & de Langre (2009), the kinetic energy of the fluid in a volume $l \times l \times H$ is

$$ \hat{\mathbf{E}}_{\text{fluid}} = \rho l^2 \int_{0}^{H} \hat{\mathbf{u}} \hat{\mathbf{u}}^* + \hat{\mathbf{w}} \hat{\mathbf{w}}^* \mathrm{d}z, \qquad (2.14) $$

where * denotes the complex conjugate of the corresponding quantity. Similarly, the perturbation energy in an individual plant is given by

$$ \hat{\mathbf{E}}_{\text{solid}} = M \hat{\xi}_x \hat{\xi}_x^* + R \hat{q}_x \hat{q}_x^*. \qquad (2.15) $$

$\hat{\mathbf{E}}_{\text{fluid}}$ and $\hat{\mathbf{E}}_{\text{solid}}$ are used to compute the fraction of the total perturbation energy stored in the oscillating canopy, which provides indications about the location of the mode:

$$ \eta = \frac{\hat{\mathbf{E}}_{\text{solid}}}{\hat{\mathbf{E}}_{\text{solid}} + \hat{\mathbf{E}}_{\text{fluid}}}. \qquad (2.16) $$

If $\eta \approx 1$ the mode is mostly a structural or canopy mode, and if $\eta \approx 0$ the mode is mostly a fluid mode.

The present model, adapted from Py et al. (2006) to account for eddy viscosity and use a smooth mean velocity profile, allows the dynamic linear stability to be simulated for the same interaction scenarios as with the LES model. However, the LES simulations must be performed first so as to extract the mean wind and SGS turbulent viscosity profiles in order to impose them on the linear stability simulations.

# 3. Main characteristics of wind-plant interaction

## 3.1. Instantaneous behaviour

Before the average characteristics of wind flow and plant motion are presented, we find it interesting to have a qualitative look at instantaneous wind-plant interactions, as simulated by the LES model. For this purpose, three figures are described in parallel. Figure 3 shows a time sequence of wind-plant interaction in a vertical streamwise plane over a 0.90 s period in the high wind speed case (Case 4). The background colour indicates the intensity of the streamwise wind velocity component, the arrows show the wind direction and the white stems sketch canopy plants. Note that, for a better visualization, angular plant displacements represented by the white stems were multiplied by a factor of 5 in figure 3. Considering the stem located at $x = 10h$ in figure 3, figure 4 presents the 7.5 s time series of (i) the wind velocity components $u$, $v$ and $w$ at the stem top (figure 4a), (ii) the stem deflection amplitudes $q_x$ and $q_y$ (figure 4b), (iii) the stem velocity components $\zeta_x$ and $\zeta_y$ (figure 4c) and (iv) the magnitude of the different terms of the stem motion equation (2.2) (figure 4d). The dashed rectangle in figure 4 indicates the time period corresponding to the snapshot sequence shown in figure 3. The sway motion of this reference stem (displacement and velocity) during a 30 s period is also presented in figure 5 in a streamwise × spanwise axes graph, as usually reported in the forestry literature from measurement of tree motions.

At $t = 0$ s, a wind gust penetrates into the canopy around $x = 8h$, inducing a forward deflection of a group of plants (figure 3a). The drag term increases first and is opposed to the inertia and stiffness terms (figure 4d). Then, with increasing plant velocity and deflection, the inertia term changes sign and augments plant deflection while the drag term starts to decrease with wind speed. At maximum plant deflection, the stiffness
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Streamwise cross-section of instantaneous wind-plant interaction at 0.30 s intervals during a period of 0.90 s (Case 4). The background colours represent the magnitude of the streamwise wind velocity, the arrows the wind vectors and the white stems the plant displacements under wind forcing. For a better visualization, angular plant displacements were multiplied by a factor of 5.</figcaption></figure>

and gravity term reach opposite maxima, and the plant velocity is zero. Just after the gust, plants spring back ($t = 0.30$ s in figure 3) and oscillate around their axis ($t = 0.60$ and 0.90 s), before their motion is damped (as shown in figure 4b from $t = 0–4$ s) and they get hit by another gust. During this period when the plant sways after the passage of the gust, the dominant terms in the plant motion equation appear to be the stiffness, the inertia and the drag terms while the gravity and damping terms are smaller (figure 4d). The signatures of plant displacement and velocity represented
<figure><img src="image_6.png" /><figcaption>FIGURE 4. Time series of wind velocity components (a), plant deflection (b), plant velocity components (c) and magnitude of the various terms of the stem motion equation (2.2) (d) at canopy top and x = 10h in figure 3 for Case 4. The dashed rectangle indicates the time period corresponding to the snapshot sequence of figure 3.</figcaption></figure>

in figure 5 illustrate the complex motion of crop plants although their mechanical
characteristics are 'simple'. These signatures have a similar shape as those recorded by
Peltola (1996) for Scots pine trees and by James, Haritos & Ades (2006) for various
tree types. We can observed from figure 5(a) that, in Case 4, plant displacements
reach up to 0.15 m for the strongest wind gusts. Consequently, the limits of validity of
the small-displacement assumptions considered in our plant motion modelling may
be reached in Case 4 for the strongest wind gusts.

Wind gusts inducing plant swaying are characterized by large positive values of *u*
and large negative values of *w* (figure 4a). This defines the signature of sweep motions,
i.e. downward motions of those coherent structures that are typical of canopy flows.
Unlike wind velocity, time series of plant deflection (amplitude and velocity) exhibit
a dominant periodicity of about 1 s (figure 4b, c), which corresponds to the natural
vibration frequency of alfalfa plant, *f*₀ = 1.05 Hz. The amplitude of plant deflection
<figure><img src="image_7.png" /><figcaption>FIGURE 5. Displacement (a) and velocity (b) components of the stem located at x = 10h in figure 3 during a period of 30 s, for Case 4.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 6. Simulated vertical profiles of mean horizontal wind velocity (a), momentum flux (b), standard deviations of the three wind components (c) ($\sigma_u$: solid line; $\sigma_v$: long dashed line; $\sigma_w$: small dashed line), skewnesses of u and w (d) ($Sk_u$: solid line; $Sk_w$: dashed line), SGS turbulent Reynolds number used by the linear stability analysis in §7 (e), for Case 3. All variables are normalized by the mean streamwise wind velocity at tree top, $U_h$, or the friction velocity above the canopy, $u_*$.</figcaption></figure>

is about 0.1 m at canopy top for a wind speed of about 10 m s⁻¹ (figure 4a, b). This value is in agreement with the 0.1 m deflection observed by Py et al. (2006) for alfalfa stems under a vertically averaged windload of 3 m s⁻¹, which corresponds, from the average wind profile presented hereafter (figure 6a), to a wind speed of about 9 m s⁻¹ at canopy top.

## 3.2. Mean flow and plant motion statistics

The basic normalized profiles of turbulent wind statistics (i.e. wind velocity $\langle u \rangle_{xyt}$, momentum flux $-\langle u'w' \rangle_{xyt}$, standard deviations of the three wind velocity components $\sigma_u$, $\sigma_v$ and $\sigma_w$, streamwise and vertical velocity skewnesses $Sk_u$ and $Sk_w$) are presented in figure 6 for Case 3 (the normalized profiles for the other cases are similar, and the corresponding figures are not shown). The subtle oscillations appearing just above the canopy on some profiles result from small numerical perturbations induced by the sharp transition between grid cells with and without vegetation; they should not impact the main flow dynamics. For the same case, figure 7 presents the distributions of plant deflection amplitude $q_i$ and velocity $\zeta_i$ in the streamwise and spanwise
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Distribution of plant deflection (a) and plant velocity (b) components at canopy top for Case 3.</figcaption></figure>

directions. Table 2 summarizes the main statistics of wind flow and plant motion at canopy top for all cases.

Turbulent wind statistics profiles exhibit the well-known behaviour of vegetated canopy flow (see Finnigan 2000 for a review), namely (i) a strong shear at canopy top associated with an inflection point in mean horizontal velocity (figure 6a), (ii) a rapid decrease within the canopy of the three wind velocity standard deviations and momentum flux (figures 6b and 6c) and (iii) positive and negative skewnesses of the streamwise and vertical wind velocity components at canopy top, respectively. This behaviour reflects that turbulence is dominated by intermittent, energetic downward-moving gusts. The shear length scale $L_s = U_h/(\partial \langle u \rangle_{xyt} / \partial z)_{z=h}$ (where $U_h$ is the mean wind velocity at canopy top), which characterizes the vertical scale of coherent structures (Raupach et al. 1996), is about 0.3h and does not depend on the wind speed at canopy top (table 2).

As one might expect, both distributions of $q_i$ and $\zeta_i$ are symmetric around their mean value in the spanwise direction (figure 7). In the streamwise direction, the distribution of $q_x$ exhibits a longer tail in forward displacements (positive values) than in backward displacements (negative values) as the latter ones are against the mean wind flow. This feature appears to be less pronounced with increasing $U_h$ as the large positive value of $Sk_{q_x}$ decreases (table 2). Furthermore, backward motions can be faster than forward motions due to the plant stiffness. This is illustrated by the negative value of the skewness of $\zeta_x$ (see table 2), indicating a slight asymmetry in the distribution of $\zeta_x$ with higher negative values than positive values. In the same way as $Sk_{q_x}$, $Sk_{\zeta_x}$ decreases with increasing $U_h$. The ratio between the standard deviations of wind and plant velocities at canopy top decreases with increasing wind speed (table 2), $\sigma_u/\sigma_{\zeta_x}$ goes from 110 to 16 and $\sigma_v/\sigma_{\zeta_y}$ from 124 to 19 with $U_h$ increasing from 1.0 to 3.8 m s$^{-1}$, indicating that wind turbulence becomes more effective in inducing plant motion at larger wind speed. For the range of wind speed considered in this study, the standard deviation of plant velocity appears still much lower than that of the wind velocity.

We saw in the last subsection that a variation in the amplitude of plant deflection with wind speed at canopy top is in agreement with the values observed in Py et al. (2006). The deflection velocity is lower than that of the wind flow by slightly more than one order of magnitude (figure 4a, c). Figure 8 shows a comparison of simulated and measured standard deviations of plant velocity $\sigma_\zeta$ (where $\zeta$ is the scalar plant velocity at canopy top), over a range of wind speed. The two sets of values are in very
<table><thead><tr><th>Variable</th><th>Designation</th><th>Case 1</th><th>Case 2</th><th>Case 3</th><th>Case 4</th><th>Case 4 bis*</th></tr></thead><tbody><tr><td><strong>Wind</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>U<sub>h</sub></td><td>Wind speed (m s<sup>-1</sup>)</td><td>1.0</td><td>2.0</td><td>2.9</td><td>3.8</td><td>3.9</td></tr><tr><td>L<sub>s</sub>/h</td><td>Shear length scale</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td></tr><tr><td>σ<sub>u</sub></td><td>Standard deviation of u (m s<sup>-1</sup>)</td><td>0.65</td><td>1.26</td><td>1.90</td><td>2.48</td><td>2.54</td></tr><tr><td>σ<sub>v</sub></td><td>Standard deviation of v (m s<sup>-1</sup>)</td><td>0.42</td><td>0.82</td><td>1.23</td><td>1.61</td><td>1.63</td></tr><tr><td>σ<sub>w</sub></td><td>Standard deviation of w (m s<sup>-1</sup>)</td><td>0.33</td><td>0.64</td><td>0.95</td><td>1.25</td><td>1.26</td></tr><tr><td>Sk<sub>u</sub></td><td>Skewness of u</td><td>0.64</td><td>0.66</td><td>0.72</td><td>0.72</td><td>0.70</td></tr><tr><td>Sk<sub>v</sub></td><td>Skewness of v</td><td>-0.04</td><td>0.08</td><td>-0.02</td><td>-0.09</td><td>-0.00</td></tr><tr><td>Sk<sub>w</sub></td><td>Skewness of w</td><td>-0.50</td><td>-0.50</td><td>-0.50</td><td>-0.50</td><td>-0.50</td></tr><tr><td>U<sub>c</sub>/U<sub>h</sub></td><td>Normalized convection velocity</td><td>1.6</td><td>1.5</td><td>1.5</td><td>1.6</td><td>1.6</td></tr><tr><td><strong>Plant</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>σ<sub>qx</sub></td><td>Standard deviation of q<sub>x</sub> (m)</td><td>0.0015</td><td>0.0066</td><td>0.0160</td><td>0.0275</td><td>0.0271</td></tr><tr><td>σ<sub>qy</sub></td><td>Standard deviation of q<sub>y</sub> (m)</td><td>0.0007</td><td>0.0034</td><td>0.0082</td><td>0.0139</td><td>0.0137</td></tr><tr><td>Sk<sub>qx</sub></td><td>Skewness of q<sub>x</sub></td><td>2.06</td><td>1.28</td><td>0.98</td><td>0.78</td><td>0.76</td></tr><tr><td>Sk<sub>qy</sub></td><td>Skewness of q<sub>y</sub></td><td>-0.02</td><td>0.20</td><td>0.03</td><td>-0.14</td><td>-0.02</td></tr><tr><td>σ<sub>ξx</sub></td><td>Standard deviation of ξ<sub>x</sub> (m s<sup>-1</sup>)</td><td>0.0059</td><td>0.0338</td><td>0.0880</td><td>0.1585</td><td>0.1543</td></tr><tr><td>σ<sub>ξy</sub></td><td>Standard deviation of ξ<sub>y</sub> (m s<sup>-1</sup>)</td><td>0.0034</td><td>0.0190</td><td>0.0482</td><td>0.0849</td><td>0.0827</td></tr><tr><td>Sk<sub>ξx</sub></td><td>Skewness of ξ<sub>x</sub></td><td>-0.43</td><td>-0.30</td><td>-0.24</td><td>-0.16</td><td>-0.16</td></tr><tr><td>Sk<sub>ξy</sub></td><td>Skewness of ξ<sub>y</sub></td><td>-0.08</td><td>-0.03</td><td>0.01</td><td>0.03</td><td>0.01</td></tr><tr><td><strong>Wind–plant interaction</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>σ<sub>u</sub>/σ<sub>ξx</sub></td><td>Standard deviation ratio between u and ξ<sub>x</sub></td><td>110</td><td>37</td><td>22</td><td>16</td><td>16</td></tr><tr><td>σ<sub>v</sub>/σ<sub>ξy</sub></td><td>Standard deviation ratio between v and ξ<sub>y</sub></td><td>124</td><td>43</td><td>26</td><td>19</td><td>20</td></tr><tr><td>R<sub>x</sub></td><td>Normalized streamwise r.m.s. difference between the drag terms at canopy top for waving and fixed plants</td><td>0.09</td><td>0.17</td><td>0.26</td><td>0.51</td><td>—</td></tr><tr><td>R<sub>y</sub></td><td>Normalized spanwise r.m.s. difference between the drag terms at canopy top for waving and fixed plants</td><td>0.18</td><td>0.44</td><td>0.59</td><td>1.03</td><td>—</td></tr></tbody></table>

<caption>TABLE 2. Main statistics of simulated wind flow and plant motion at canopy top.</caption>

good agreement. They both increase with $U_h$ due to the enhancement of turbulence induced by the larger wind shear at canopy top (see the values of the three standard variations of wind velocity components in table 2). Although plant deflections are small we need to remember here that the assumption of linear deformation of the plant stem considered in our model may induce a slight underestimation of simulated canopy-top plant deflection and velocity compared to a flexible plant, and this underestimation should increase with wind speed. Regarding the impact of not accounting for streamlining effect in our model, it is difficult here to evaluate its consequences since, on one hand, the velocity of penetrating wind gusts within the canopy should increase with plant deformation while, on the other hand, the plant frontal area density should be reduced, and so the drag force.

In conclusion, the visualization of instantaneous wind–plant interactions simulated by our model over an alfalfa crop canopy confirms the realism of the model despite the simplifications considered in our model. Although the limits of validity of the small-deflection assumption used in our model may be reached for some strong wind gusts in Case 4, the variations in the magnitude and velocity of plant displacements

<footnote>* Same as Case 4 but with fixed plants instead of waving plants.</footnote>
<figure><img src="image_10.png" /></figure>

FIGURE 8. Comparison between measured (empty circles) and simulated (black squares) standard deviation of alfalfa plant velocity against the wind speed at canopy top. The experimental dataset comes from Py *et al.* (2006).

<figure><img src="image_11.png" /></figure>

FIGURE 9. Snapshot of the simulated alfalfa crop motion for Case 4. For a better visualization, angular plant displacement was multiplied by a factor of 5.

with wind speed at canopy top are in very good agreement with *in situ* measurements performed by Py *et al.* (2006) over a similar crop canopy.

4. Plant waving structures

After it was verified that the LES model accurately simulates plant deflection and velocity over a current range of wind velocity, we now focus on the main characteristics of plant waving structures. Figure 9 shows an instantaneous three-dimensional view of the simulated crop motion. Plant displacements have been accentuated in order to have a better view of the waving structures. The dark patches appearing on the crop surface correspond to regions where plants are strongly deflected downwind under
the action of strong sweeping wind gusts, as those observed in the previous section from the sequence of wind-plant interactions in a vertical plane orientated streamwise (figure 3). It can be seen by looking at animations of crop motion that these patches move essentially along the mean wind direction. They are induced by the development or impingement of coherent structures at canopy top. The relationship between crop waving structures and coherent eddies is investigated in § 5. As was stated in § 1, these wave-like crop motions are known as ‘honami’ waves, and resemble the cat’s paws patterns observed on water surfaces. These patches correspond to the white patterns observed by Py et al. (2005) on their alfalfa field (figure 1a).

In the same way as was done by Py et al. (2005, 2006), the main spatio-temporal features of these plant waving structures were extracted in all four cases from a BOD of the crop velocity fields $\zeta(x, y, t) = [\zeta_x, \zeta_y]$, which were recorded at 10 Hz during a 30 s period. The BOD approach was first introduced by Aubry, Guyonnet & Lima (1991). The reader can refer to Hémon & Santi (2003) for a complete review of the approach and to Py et al. (2006) for its application to plant motion. To summarize, the BOD allows $\zeta(x, y, t)$ to be decomposed into a finite series of spatio-temporal structures as follows:

$$ \zeta(x, y, t) = \sum_{k=1}^{N} \sqrt{\alpha_k} \mu_k(t) \psi_k(x, y), \quad (4.1) $$

where $\mu_k$ and $\psi_k = [\psi_{kx}, \psi_{ky}]$ are, respectively, the temporal and spatial functions of mode $k$, referred to as ‘chronos’ and ‘topos’ and $\sqrt{\alpha_k}$ is the weight factor of each spatio-temporal structure ($\mu_k, \psi_k$). ‘Chronos’ and ‘topos’ form a set of orthogonal functions, that are the eigenfunctions of the temporal and spatial correlation operators of $\zeta_i$ with the same eigenvalues $\alpha_k$, respectively. As explained by Py et al. (2005), compared to other decomposition approaches such as the empirical orthogonal functions (EOF) or the proper orthogonal decomposition (POD), the BOD has the advantages of being applicable to space-time signal without any assumption other than being square-integrable, while POD and EOF approaches require also the ergodicity, the stationarity and a Gaussian distribution of the signal, which makes them not applicable in heterogeneous conditions or for signals with intermittent events. But the most important difference for our study is the fact that BOD performs an analysis of the signal in both space and time, allowing to extract spatial and time information on main coherent structures of crop-plant motions.

The spatio-temporal modes ($\mu_k, \psi_k$) are ranked according to the descending order of their kinetic energy $\alpha_k$. For Case 1 (but similar results are observed for the other cases), figure 10 presents the cumulative energy, as recovered by the BOD, of the crop velocity field as a function of the rank of the spatio-temporal modes. Similarly to the crop motion analysis of Py et al. (2005) we observe a good convergence of the BOD, as 75% of the energy is reproduced in the signal by the first 20 spatio-temporal modes. This rapid convergence of the BOD indicates the presence of large coherent structures in the crop velocity field, whose temporal and spatial characteristics are defined by the first sets of ‘chronos’ and ‘topos’, respectively. Hence we only focus hereafter on the most two energetic modes.

The divergence of the spatial eigenvectors ($\partial\psi_{kx}/\partial x + \partial\psi_{ky}/\partial y$), i.e. the ‘topos’, of the most two energetic modes and their associated temporal eigenvector, i.e. the ‘chronos’, are presented in figures 11 and 12 for the lowest (Case 1) and highest (Case 4) wind speed. The amplitudes of both spatial and temporal eigenvectors are not shown since only their patterns are important here. In both cases organized motions clearly appear
<figure><img src="image_12.png" /><figcaption>FIGURE 10. Percentage of the total kinetic energy recovered versus the number of spatio-temporal modes considered in the BOD of the plant velocity field in Case 1.</figcaption></figure>

from ‘topos’ as large parallel stripes perpendicular to the mean wind direction. The wavelength $\lambda_p$ of these stripes is smaller in Case 1 than in Case 4 (figures 11a and 12a), i.e. $\lambda_p$ increases with mean wind speed. The wavelength $\lambda_p$ was deduced from the averaged two-dimensional Fourier transform of $\psi_k$ ($k=1, 2$) (figures 11c and 12c). While a unique wavelength peak is observed in Case 1 around 2.0h (figure 11c), two distinct peaks of similar magnitude appear in Case 4 (figure 12c), one around 4.4h and the other one around 8.7h. This feature can be explained by the presence of high wind velocity regions (about 5.0 m s$^{-1}$) between $y=0$ and 10h, elongated in the streamwise direction, and lower velocity regions (about 2.5 m s$^{-1}$) between $y=15$ and 20h. As it will be discussed in the next section, such meandering elongated structures are typical of logarithmic region of neutral boundary layers.

The lifetime of these elongated structures is larger than the application time of the BOD (30 s). Consequently, the 4.4h wavelength peak is associated with the low wind speed region and the 8.7h with the high wind speed region. A smaller wavelength is indeed perceptible in figure 12(a) between $y=15$ and 20h. Such structures are also present in other simulated cases but the difference in wind speed between the high and low wind speed regions are too small for two distinct wavelength peaks to be observed in the Fourier transform. These spatial organized motions are associated with a regular oscillating behaviour appearing on ‘chronos’. The averaged Fourier transform of $\mu_k$ ($k=1, 2$) (figures 11d and 12d) indicates a well-defined common frequency $f_p \approx 1.05$ Hz in Cases 1 and 4, as in the other cases too (not shown), as was previously observed from plant velocity time series (figure 4c). This frequency is equal to the natural vibration frequency $f_0$ of alfalfa plants. Both ‘topos’ and ‘chronos’ of the first and second modes appears phase-lagged in space and in time, respectively. As a consequence, these organized motions propagate along the main wind direction with a phase velocity $U_{cp} = \lambda_p f_p$.

The normalized wavelength $\lambda_p/h$, frequency $f_p/f_0$ and phase velocity $U_{cp}/U_h$ of the main spatial organized structures deduced from ‘topos’ and ‘chronos’ are plotted in figure 13 against the reduced velocity $U_r = U_h/(f_0h)$, $f_0h$ being a characteristic velocity of plant stem. Compared with values deduced by Py *et al.* (2006) from video-recorded alfalfa crop motion, the simulated values of $\lambda_p/h$, $f_p/f_0$ and $U_{cp}/U_h$ are in fairly
<figure><img src="image_13.png" /><figcaption>FIGURE 11. First and second 'topos' (spatial eigenvector) divergence and 'chronos' (temporal eigenvector) of alfalfa plant velocity field in Case 1 (a and b, respectively). Average 'topos' (c) and 'chronos' (d) spectra of the first two modes are shown on the right-hand side.</figcaption></figure>

good agreement although $\lambda_p/h$ and, consequently, $U_{cp}/U_h$, are slightly overestimated by the model. The ratio $\lambda_p/h$ appears to increase with $U_r$ while $f_p/f_0$ is independent of $U_r$. The phase velocity $U_{cp}$ of crop motion is around $1.4U_h$, except in the lower wind speed case where it is slightly larger ($2.0U_h$). These values of $U_{cp}$ are consistent with the average value of $1.6U_h$ observed by Finnigan (1979) over a uniform wheat canopy. The reason for the slight overestimation of $\lambda_p/h$ by our model is not clear but it may be related to the homogeneity of alfalfa plant properties considered in our simulation, as compared with the variability of plant properties in a real crop. This variability inside the crop concerns in particular the natural vibration frequency $f_0$, the height $h$ and the mass $m$ of plants, for which Py et al. (2006) observed a range of values of about 0.8–1.5 Hz, 0.47–0.84 m and 0.0039–0.0186 kg, respectively. The simplifications considered in our plant model should have a small impact on the main characteristics of coherent waving structures since $\lambda_p$ depends mostly on $f_0$ and on the convection velocity of canopy-top wind gusts, and less on plant deflections, as
<figure><img src="image_14.png" /><figcaption>FIGURE 12. Same as figure 11 but for Case 4.</figcaption></figure>

it will be seen in § 6. The convection velocity of coherent structures may depend on plant deformation as it should increase with more streamlined plant shape. However, as already mentioned, this streamlining effect should be limited in our cases due to the small values of plant deflections.

To conclude, our model appears to simulate accurately the main coherent motions of the crop canopy, as compared with the video recordings by Py *et al.* (2006). These coherent motions are characterized by a frequency close to the natural vibration frequency of the plants and by a wavelength that increases with wind speed at canopy top.

**5. Coherent crop motion and coherent eddy structures**

The organized crop motion or ‘honami’ waves identified in the previous section are initiated by the development or impingement of coherent eddy structures at canopy-top, but the nature of their interaction is still largely unknown. In the present section we investigate the differences between the main properties of canopy-top coherent eddy structures and coherent crop motion. But, first, we find it interesting to look
<figure><img src="image_15.png" /><figcaption>FIGURE 13. Comparison between experimental observations (empty circles), LES (black squares) and linear stability analysis (empty triangles) of the normalized wavelength (a) frequency (b) and phase velocity (c) of coherent waving patches (extracted from BOD of alfalfa plant velocity field for measurements and LES) versus the reduced velocity. The experimental dataset comes from Py et al. (2006).</figcaption></figure>

at instantaneous flow fields above the canopy in order to visualize large outer-layer structures compared to canopy structures.

5.1. *Outer-layer structures above the canopy*

Figure 14 shows contours of streamwise fluctuations at a given time, above ($z=3h$) and at canopy top, for Case 4. Streamwise fluctuations at canopy top are characterized by small longitudinal patterns that may be the signature of coherent structures induced by the canopy itself. These structures should be related to plant motions patterns identified in figure 9. With increasing height, these structures increase in size and look like elongated streamwise structures. Such structures have been previously observed experimentally and numerically, for neutral stratification, in logarithmic region of near-wall boundary layers and in the atmospheric surface layer (see e.g. Moeng & Sullivan 1994; Kim & Adrian 1999; Drobinski et al. 2004; Foster et al. 2006; Hutchins & Marusic 2007). In the atmosphere, these elongated structures are known as steak structures (Drobinski et al. 2004). Their size is usually related to the boundary layer thickness $\delta$ (Hutchins & Marusic 2007) as well as to the surface
<figure><img src="image_16.png" /><figcaption>FIGURE 14. Snapshot of horizontal cross-sections (x - y) of the streamwise wind velocity fluctuations at z = 3h (a) and z = h (b), for Case 4.</figcaption></figure>

roughness length (Lin *et al.* 1997). They can extend to over 20δ in length (Hutchins & Marusic 2007), and their width was observed by Lin *et al.* (1997) and Drobinski *et al.* (2007), from LES of the planetary boundary layer, to increase with height and to be of the order of the atmospheric surface layer thickness. Streaks present in the atmospheric surface layer have an average spacing of hundred of meters and a time scale of several minutes (Drobinski *et al.* 2004). Their spatial size and lifetime are therefore much larger than our crop size and the duration of our simulations, respectively. Consequently, we are not expected in our domain to simulate near-surface streaks, as already mentioned in §2.3. Such structures should only be seen here as background wind flow. The elongated structures observed in figure 14 are locally induced by the wind shear in the logarithmic region located above the canopy. These structures may be considered as part of streak structures, which could be consistent with the suggestion of Adrian, Meinhart & Tomkins (2000) and Hommema & Adrian (2003) that streaks may be associated with packets of hairpin vortices.

In our simulations, longitudinal structures reach about 4h in width at canopy top and 8h in width at z = 3h. Removing the depth of the Rayleigh damping layer (2.5 m), the simulated boundary layer has a depth of about 5.5 m (i.e. 8h). Hence, structure widths correspond to 0.5δ and 1δ, at canopy top and z = 3h, respectively, which are in the range of expected values. The limited size of our computational domain and
the use of periodic conditions may impact the realism of large structures simulated in the upper layers of the domain. However, since the turbulent structures of interest in our study are canopy structures that scale with the canopy height, the size of our domain should be sufficient for accurately simulating such structures as previously demonstrated from LES with similar domain size compared to canopy height (see e.g. Su et al. 1998, 2000; Watanabe 2004; Dupont & Brunet 2008a, 2009). The interaction between canopy structures and large outer-layer structures is still a research issue. Hunt & Morrison (2000), Hunt & Carlotti (2001) and Carlotti (2002) suggested that a 'top-down' trajectory mechanism for main coherent structures within the atmospheric surface layer may occur. In other words, outer-layer eddies may impinge onto the ground, inducing the development of internal boundary layers, as they are blocked by the presence of the ground, in which small eddies or surface eddies developed (Fesquet et al. 2009). In our case over an homogeneous canopy, structures located in the roughness sublayer, say between $z=h$ and $2h$, may impinge onto the canopy top and initiate the development of canopy structures while large outer-scaled structures may only be seen as a footprint on the canopy-top mean wind flow and plant motion, as observed in the last section on 'topos' of plant motion for Case 4.

## 5.2. Canopy structures

Coherent eddy structures over homogeneous vegetation canopies have been investigated for years from *in situ* and wind-tunnel experiments as well as from LES. They scale with mean canopy height and have a convection velocity $U_c$ of about $1.8U_h$ (Shaw et al. 1995; Finnigan 2000). The development, characteristics and length scales of these coherent structures appear similar to those observed in plane mixing-layer flows (Raupach et al. 1996). In the light of this analogy Raupach et al. (1996) deduced from many field and wind-tunnel datasets that the mean longitudinal separation or wavelength $\lambda_w$ between adjacent coherent structures is only function of the shear length scale $L_s$ (defined in §3.2) and independent of wind speed, as for plane mixing-layer flows.

Using the spatial velocity correlations deduced from the simulations, we verified that the spatial scales of coherent eddies at canopy top are in agreement with previous observations (Shaw et al. 1995; Su et al. 2000; Dupont & Brunet 2008a). The average convection velocity $U_c$ of coherent eddies at canopy top, deduced from space-time correlation of the vertical wind velocity $w$ is about $1.55U_h$ (table 2), which is in relatively good agreement with the usual value of $1.8U_h$ considered in the literature. This value is also close to that of the phase velocity of coherent crop motion, $U_{cp} = 1.4U_h$, observed in §4.

In order to identify the wavelength $\lambda_w$ of wind-flow coherent structures, we first attempted to perform a BOD of $w$, as was done for the crop velocity field in the previous section. The vertical wind velocity component $w$ was preferred to the streamwise component $u$, as the latter includes contribution from large eddies coming from the atmosphere above, thereby making $w$ more representative of the active turbulence at canopy top (Raupach et al. 1996). However, the BOD of $w$ did not allow dominant energetic modes to be identified. More than 100 modes were necessary to recover 75% of the signal kinetic energy, this feature being probably due to the large background turbulence surrounding coherent structures. For this canopy wind flow, the BOD degenerates into a Fourier decomposition as observed in homogeneous turbulence (Farge et al. 2003). For this reason, $\lambda_w$ was instead deduced in the four simulated cases from the peak wavenumber of the resolved-scaled spectrum of $w$ at canopy top, averaged over a 30 s period (figure 15).
<figure><img src="image_17.png" /><figcaption>FIGURE 15. Average resolved-scaled spectra of the vertical wind velocity at canopy-top for the four simulated cases (Case 1: dashed-dot line; Case 2: dashed line; Case 3: dashed-dot-dot line; Case 4: solid line).</figcaption></figure>

The overall magnitude of w-spectra increases with increasing wind speed due to the enhancement of TKE at canopy top (figure 15). In all cases the w-spectra display the same familiar shape of canopy top w-spectra with a $k_1^{+1}$ slope at low wavenumbers and a slightly steeper slope than the $-2/3$ power law spectral subrange at intermediate wavenumbers. This shape of w-spectra is consistent with, for example, measurements from Amiro (1990) over a spruce canopy or LES results from Su et al. (1998) over a forest canopy. This steeper spectral slope observed within canopies compared to the atmospheric surface layer, is usually viewed as an evidence of short-circuit of the inertial cascade due to the interaction of large eddies with vegetation elements (Finnigan 2000). It may also be induced by a too dissipative subgrid scale model within the canopy (Su et al. 1998) At high wavenumbers the spectrum decays at a much larger rate than $k_1^{-2/3}$, due to the limited spatial resolution of the domain and to the SGS model that are felt at smaller wavenumbers than the cutoff wavenumber. This discrepancy is typical of LES (Carlotti 2002; Foster et al. 2006). In the four cases, the spectra exhibit the same peak at around $4h \pm 0.5h$. Kaimal & Finnigan (1994) reported that the spectral peak frequency $f_w$ of w measured at canopy top has been observed as identical over a large range of canopies, including crops and forests, and equal to $f_w h/U_h \approx 0.45 \pm 0.05$. As in our case the average convection velocity $U_c$ of turbulent structures is about $1.55U_h$, the spatial wavelength $\lambda_w$ is around $3.5h \pm 0.4h$ (using Taylor's frozen turbulence hypothesis). Thereby, the position of the w-spectral peak around $4h$ and its independence from wind speed is in excellent agreement with measurements reported by Kaimal & Finnigan (1994). However, $\lambda_w$ appears larger than the value of $2.4h$ predicted by the plane-mixing layer analogy of Raupach et al. (1996) ($8.1L_s$ with $L_s = 0.30h$ here).

Even though the wavelengths of organized crop motion $\lambda_p$ and turbulent structures $\lambda_w$ are of the same order, they behave differently with increasing wind speed $U_h$. The wavelength of crop motion increases with $U_h$ while its frequency is independent of $U_h$. Conversely, the wavelength of the turbulent coherent structures does not depend on $U_h$, as is predicted by the plane-mixing layer analogy, but their frequency increases with $U_h$. Consequently, organized crop motions cannot be considered as
direct signatures of coherent eddy structures of the flow, although the latter may initiate the former. Wavelengths of both structures are close to each other for a reduced velocity of about 4, which means that alfalfa plants may be in phase with the coherent structures of the flow for a mean wind speed of about 2.9 m s⁻¹ at canopy top. A resonant interaction is therefore possible but probably attenuated by the background turbulence of the flow.

The possibility of deducing the main characteristics of the wind flow from the velocity field of the crop motion (e.g. from video recording) does not appear straightforward. Nevertheless, it can be imagined that in the future (i) the knowledge of the relationship between $λ_p$ and $U_r$ as well as the knowledge of plant mechanical properties may give access to the mean wind speed at canopy top, and (ii) the TKE of the flow at canopy top may be deduced from the plant velocity variance $\sigma_\zeta$, as $\sigma_\zeta$ was observed to increase with $\sigma_u$. Reconstructing the time and spatial wind field at canopy top from the crop velocity field, by using for example a crop motion model in an ‘inversion’ mode, remains challenging and is out of the scope of the present paper.

## 6. Impact of crop motion on turbulent wind flow

In order to estimate the possible impact of plant motion on the wind flow, we computed the non-dimensional spatial average of the root mean square (r.m.s.) value of the difference between the canopy drag in direction $i$ of a moving canopy and that of a fixed one, $R_i$ ($R_1 = R_x$, $R_2 = R_y$), as follows:

$$R_i = \left\langle \left[ \left| u_i - z \frac{\partial q_i}{\partial t} \right|^2 \left( u_i - z \frac{\partial q_i}{\partial t} \right)^2 - |u_i|^2 u_i^2 \right]^{1/2} / \left| u_i \right| \right\rangle_{xyt}^{z=h}, \quad (6.1)$$

where $i \in \{1, 2\}$. A zero value of $R_i$ means that plant motion has no impact on the wind flow. Values of $R_i$ obtained in the four simulated cases are given in table 2. As for the variance of the crop velocity field, $R_i$ increases with increasing wind speed, meaning that the impact of plant motion on wind flow is enhanced, as was also deduced by Finnigan & Mulhearn (1978b) from their analytical model. However we do not know at this stage whether the effect of plant motion is negligible or not on the wind flow. For this reason one additional simulation was performed with fixed plants in conditions similar to those of Case 4, for which the $R_i$ components are the highest. This simulation is referred as Case 4bis. No differences on the canopy-top main statistics of the wind flow and plant motion (see table 2), as well as on the basic wind profiles and the coherent structure wavelength $\lambda_w$, were observed between waving and fixed plants (figures not shown). Consequently, plant motion appears to have negligible effect on the wind flow in the range of usual wind speed values considered in this study. The much lower values of the magnitude and variance of plant velocity (see table 2) compared to wind flow ones explain certainly the negligible effect of plant motion on the wind flow.

In order to further investigate the impact of plant motion on the wind flow, we performed a sensitivity study of the streamwise r.m.s. difference $R_x$ (see (6.1)) to the main crop plant properties. The goal here is to see how the modification of the main plant properties may increase the velocity of plants compared to that of the wind flow, and consequently may enhance the possible influence of plant motions on the wind flow. Due to considerations regarding computational time, we then used a simplified version of our model. This consisted in solving the crop motion equation (2.2) only for the streamwise component ($i=1$) and forcing this equation with a continuous
<table><caption>TABLE 3. Values of the plant mechanical properties used in the sensitivity study presented in figure 16. The values of the four parameters for the coded value 0 correspond to the values used in the reference simulation (table 1).</caption><thead><tr><th>Parameter</th><th>Designation</th><th>-1</th><th>-0.5</th><th>0</th><th>+0.5</th><th>+1</th></tr></thead><tbody><tr><td>f<sub>0</sub></td><td>Natural vibration frequency (Hz)</td><td>0.10</td><td>0.575</td><td>1.05</td><td>1.525</td><td>2.0</td></tr><tr><td>l</td><td>Mean plant spacing (m)</td><td>0.01</td><td>0.03</td><td>0.05</td><td>0.07</td><td>0.09</td></tr><tr><td>m</td><td>Plant mass (kg)</td><td>0.005</td><td>0.0094</td><td>0.0138</td><td>0.0182</td><td>0.0226</td></tr><tr><td>ξ</td><td>Damping coefficient</td><td>0.</td><td>0.0446</td><td>0.0875</td><td>0.1321</td><td>0.1750</td></tr></tbody></table>

series of instantaneous wind profiles. These profiles were constructed from a 30s recording of the streamwise wind velocity component simulated by ARPS at the top of the alfalfa canopy in high wind speed (Case 4). The wind velocity within the canopy was then deduced from its canopy-top value by using an exponential decrease $u(z) = u(h) \exp(LAI(z/h - 1))$. With this simplified approach, plant motion cannot influence the wind flow but the $R_x$ parameter can provide useful information on the potential impact of plant motion on the wind flow. It was verified that the main properties of the wind profile (intensity and variance) and plant motion (magnitude and variance of plant displacement and velocity) obtained from this simplified model were similar to those obtained with the complete model over the alfalfa canopy.

We consider crop canopies with identical height ($h=0.69$ m) and foliar density (LAI=3) as for the alfalfa canopy studied previously. We modify in turn the plant mass $m$, the natural vibration frequency $f_0$, the damping coefficient $\xi$ and the plant spacing $l$. The range of variation of these four parameters is given in table 3, where the coded values -1 and +1 refer to the minimum and maximum values of each parameter, and the coded value 0 refers to the reference values used in the previous simulations. The maximum values are approximately the double of reference values, and the minimum values have been chosen such as the reference values are the average values between the maximum and minimum ones. The range of values of these four parameters considered here does not cover the all range of values encounter in crop plants but this range should be sufficient to give us information on how the magnitude of plant velocity compared to that of the wind flow varies with these parameters, and so on the possible influence of plant motion on the wind flow. Figure 16 presents the variation of $R_x/R_x^0$ (where $R_x^0$ refers to $R_x$ for the coded value 0) versus the coded values of the four mechanical parameters of the plant. It follows that for a given LAI, the effect of plant waving on the wind flow is enhanced with decreasing natural vibration frequency, mass and damping coefficient and increasing plant spacing (i.e. denser plants as LAI remains constant). For the latter three variables this finding is consistent with the results obtained by Finnigan & Mulhearn (1978b) with their analytical model.

We then consider a wheat plant forced by the same wind flow as for the alfalfa plant. The mechanical properties of the wheat plant are deduced from Py et al. (2006) and given in table 1. As the height, leaf-area index and vertical foliar distribution of a wheat canopy are close to those of an alfalfa canopy, it was assumed that the wind profile is similar within both canopies. The value of $R_x/R_x^0$ obtained for the wheat plant appears lower than that of alfalfa plant (figure 16). Consequently, wheat plant motion should also have a negligible impact on the wind flow over the range of wind speed considered in this study.
<figure><img src="image_18.png" /><figcaption>FIGURE 16. Sensitivity to the plant mechanical properties of the normalized streamwise r.m.s. difference between the drag terms at canopy top for waving and fixed plants. On the horizontal axis the coded values of plant mechanical properties refer to the values indicated in table 3. Each plain dot represents one run where only one parameter (mass *m*, natural vibration frequency *f*₀, damping coefficient *ξ* and plant spacing *l*) was modified compared to others that were kept at their reference value (coded value equal to 0).</figcaption></figure>

A consequence of this independence of the turbulent flow from plant motion is that the increase of crop motion wavelength with wind speed obtained from the LES model results from passive motion: a larger group of plant is deflected by passing coherent eddies as they travel faster (higher wind speed), while the frequency of plant vibration remains unchanged. This means that the main characteristics of coherent crop motion depends mostly on natural plant vibration *f*₀ and on the convection velocity of canopy-top wind gusts. This passive behaviour of plant appears in contradiction with that from the wind-tunnel experiment of Finnigan & Mulhearn (1978a) over a flexible canopy, and with the lock-in mechanism predicted by the linear stability analysis of Py *et al*. (2006). In the first case Finnigan & Mulhearn (1978a) justified the importance of accounting for plant velocity in the momentum drag force term by the fact that the variances of the streamwise velocities of wind and plant have the same order of magnitude. However we observed in our simulations that *σ*<sub>u</sub> is more than 16 times as large as *σ*<sub>ξx</sub> (table 2), which was also verified experimentally by Py *et al*. (2006). With a similar analysis to that of Finnigan & Mulhearn (1978a), we can then conclude that the effect of plant motion on the wind flow is negligible in our case. Regarding now the absence of lock-in in the LES results, the discrepancy between the LES and the linear stability analysis is investigated in the next section.

## 7. Linear stability analysis and LES

In the previous section we deduced from the LES results that plant motions have a negligible impact on the turbulent wind flow. This challenges the lock-in mechanism predicted by Py *et al*. (2006) from their linear stability analysis over a similar crop canopy. In order to understand this discrepancy between the two approaches we performed linear simulations with the model described in §2.4, for the four cases previously studied with the LES.

Linear simulations were performed with a freely-oscillating canopy and a non-moving one. The base state of the linear model is defined by the mean wind velocity
<figure><img src="image_19.png" /></figure>

FIGURE 17. Variation with the reference velocity of the dimensionless wavelength (a), frequency (b), growth rate (c) and canopy energy fraction (d) of the most unstable mode of the flow over an alfalfa canopy, as deduced from the linear stability analysis. The simulations are performed with a freely-oscillating canopy ( — — ) and a non-moving one ( - - - ).

profiles 〈u〉xyt/U_h and the mean SGS turbulent Reynolds number profile Re_t = U_h h/v_t
(where v_t = 0.1√⟨e⟩xyt√ΔxΔyΔz) extracted from the LES and shown in figures 6(a)
and 6(e), respectively. As both profiles were found almost independent of the flow
velocity, only the profiles of Case 4 (U_h = 3.8 m s⁻¹) were used in all linear stability
simulations. The plant properties are taken identical to those used in the LES (table 1).
In the same way as in Py et al. (2006) we study the temporal stability of each flow
configuration over a range of reduced velocity U_r = U_h/f_0 h. For every value of U_r
we consider the most unstable mode, i.e. the mode defined by its wavenumber and
complex frequency that has the largest positive imaginary frequency ω_j.

Figure 17 shows the variation with reduced velocity of the dimensionless wavelength (a), frequency (b) and growth rate (c) of the most unstable mode of both the freely-oscillating (solid line) and the non-moving (dashed line) configurations. Figure 17(d) shows the fraction of energy of the mode which is concentrated in the oscillating canopy (see (2.16)). The canopy movements have a distinct effect on the most unstable mode. As observed in Py et al. (2006) and Gosselin & de Langre (2009), a lock-in phenomenon is clearly visible for the freely-oscillating configuration. In the reduced velocity range, $3 < U_r < 7$, corresponding to wind velocities of $2.2 \text{ m s}^{-1}$ < $U_h < 5.1 \text{ m s}^{-1}$, the instability wavelength deviates from that of the non-moving configuration, the frequency locks onto that of the oscillating canopy $f \sim f_0$ and the growth rate surges. However, for a large range of reduced velocities, the most unstable mode of the system is mostly a fluid mode. For small reduced velocities, before lock-in begins at $U_r = 3$, the dynamics of the wind and the canopy are only
weakly coupled. The oscillation of the canopy represents less than $\eta < 0.5\%$ of the energy of the most unstable mode (figure 17d). In the lock-in range, $\eta$ increases to a little less than 5% but falls back below 1% beyond $U_r = 7.5$. Since the most unstable mode is mostly a fluid mode, it is likely that one or more other modes containing more energy in the structure would dominate the canopy response in a nonlinear coupled system, even more so outside the lock-in range.

Outside the lock-in range, as the most unstable mode is essentially a fluid mode, the main characteristics of the instability predicted by the linear stability analysis should be compared with those of coherent eddy structures simulated by the LES model. The linear theory predicts a most unstable mode of wavelength $\lambda = 4.2h$, that is constant outside the lock-in range. This agrees with the wavelength $\lambda_w = 4h$ of coherent eddy structures being independent of the flow velocity (see the results from the LES model in § 5). However, the phase velocity of the instability predicted by the linear stability model appears smaller than the convection velocity of coherent eddy structures deduced from LES, i.e. $1.00U_h$ against $1.55U_h$, respectively.

Inside the lock-in range the main characteristics of the instabilities should be representative of plant motion as the ‘fluid’ instability locks onto plant motion. The wavelength, frequency and phase velocity of the most unstable mode predicted with the linear stability analysis are compared in figure 13 with those extracted with the BOD of the alfalfa canopy motions, either video-recorded in the experiments of Py et al. (2006) or simulated with the LES model. The linear analysis predicts a wavelength which increases almost linearly with the reduced velocity. The LES also does so but the predicted values are significantly different. The frequency predicted by the linear theory matches the LES and the experiments, as all three predict motions occurring at the natural frequency of the plants. As measured in the experiments, the linear theory predicts a phase velocity around the value of the wind velocity at canopy top. As mentioned above, outside the lock-in range the behaviour of the most unstable mode of the linear analysis is quite different from the LES and the experiments. This is due to the fact that in the experiments and the LES, the wavelength, frequency and velocity are extracted from BOD performed on the canopy movements, while in the linear theory the most unstable mode is almost a pure fluid mode.

Hence, according to the linear theory, the LES Cases 3 and 4 with mean flow velocities of $2.9 \text{ m s}^{-1}$ and $3.8 \text{ m s}^{-1}$ should clearly be in the lock-in range, but their most energetic wavelengths are unaffected by crop motions, as seen in the previous section. The reasons for this discrepancy between both approaches are somewhat unclear at this stage, only possible guesses are given here. First, a similar representation of crop plant as rigid stem, with the same mechanical properties, was used in both approaches. Consequently, the simplifications that have been considered in the LES approach cannot be the reason for the absence of lock-in mechanism. The only difference in plant representation between both approaches is the two-dimensional horizontal plant motions considered in the LES compared to the one-dimensional plant motion considered in the linear approach. This latter one-dimensional motion could emphasize the interaction between wind gusts and longitudinal plant motions. Second, the linear stability analysis estimates the instability of perturbations applied on a base state flow. These perturbations should be smaller than the base state values. However, in wind-flow and plant-motion dynamics, the perturbations have been observed in LES and experiments to be of the same order as the mean values (see e.g. in table 2 the values of the standard deviations of the canopy-top wind velocity compared to its mean value). This last point could indicate a limitation of applicability of this approach. Overall, the main possible explanation for this discrepancy between
LES and linear theory may be related to the presence of a nonlinear saturation mechanism in the LES, that favours always the same wavelength independently of canopy motion. In order to verify this last explanation, a future study could consist in (i) studying with the LES the transient growth of the Kelvin-Helmholtz instability and its saturation over a canopy and (ii) developing a weakly nonlinear stability analysis in order to confirm the importance of nonlinear terms in minimizing the lock-in mechanism.

## 8. Conclusions

A novel three-dimensional model that includes a two-way coupling between the wind flow and crop plant motion has been presented here. This numerical model consists of an atmospheric LES model coupled with a simple mechanical oscillator equation for crop plant motion. The canopy is represented as a poroelastic continuous medium, which is similar in its discrete form to an infinite row of identical oscillating stems. For this initial version of the coupled model, we only considered one linear mode of plant vibration, small displacements of plants and consequently no interaction between neighbouring plants and no streamlining effect.

This model has been validated successfully against video-recorded measurements previously performed by Py et al. (2006). The magnitude of plant displacement and velocity were simulated accurately by our model for various wind speed values. The BOD of the plant velocity field has revealed the presence of coherent waving motion, whose spatial and temporal characteristics agree with those obtained by Py et al. (2005). These structures correspond to the ‘honami’ motion usually observed over cereal crops on windy days. Their frequency is close to the natural vibration frequency of the plants and their spatial wavelength is about a few canopy heights and increases with wind speed. Although these organized structures of crop movements are initiated by coherent eddy structures in the wind flow, their spatial and temporal characteristics differ. Indeed, the spatial wavelength of coherent eddy structures at canopy top were observed independent of wind speed, as predicted by the plane-mixing layer analogy of canopy flow (Raupach et al. 1996), while their frequency increases with wind speed. Hence, organized crop movements are not mere footprints of coherent eddy structures. Extracting wind-flow information from crop motion data is therefore not straightforward.

Over the range of wind speed considered in this study (1–4.0 m s⁻¹), which corresponds to typical values over crop fields, we observed that alfalfa plant motion has a negligible impact on the wind flow and that no lock-in mechanism occurs between the wind flow and plant motion. We further deduced from a simplified version of the model that similar independence of the wind flow from plant swaying should be observed for a wheat crop canopy under current wind speed. This result demonstrates that it is reasonable to consider fixed plants in wind-flow modelling at canopy scale, as is usually done.

However, an improved version of the linear stability analysis of Py et al. (2006) and Gosselin & de Langre (2009), including a realistic wind profile and an eddy viscosity deduced from LES, still turns out to predict a lock-in phenomenon in the velocity range where coherent canopy motions are observed. This discrepancy between linear stability analysis and LES may be attributed to the presence of a nonlinear saturation mechanism in LES, independent on canopy motion, which is not considered in the linear stability analysis. The next steps of this study should consist in (i) studying with the LES the transient growth of the Kelvin-Helmholtz instability and its saturation
over a canopy and (ii) developing a weakly nonlinear stability analysis in order to confirm the importance of nonlinear terms in minimizing the lock-in mechanism.

Despite the simplification considered in our plant crop model, we believe that the general picture of the wind-plant interaction that comes out from this study would still be valid with a more detailed plant representation. This is mostly motivated by the fact that these simplifications are mostly valid for small plant displacements as it was considered in the present study. Considering flexible bending stems that simply flex at their base instead of flexible stems that bend throughout their height, may slightly underestimate the magnitude of head plant displacement and velocity but the general dynamics of individual plants as well as of the crop should be well reproduced with this simplification. Finally, the discrepancy on the lock-in mechanism between the LES and linear stability analysis cannot be attributed to these simplifications since both approaches used the similar representation of crop plants.

The LES model can be helpful to better understand plant canopy motion, at any position in a given landscape. It can provide useful information on (i) the impact of plant motion on plant growth (thigmomorphogenesis) and (ii) plant vulnerability to windload. The LES model was only applied here to crop canopies with a simple plant structure, as compared with that of a tree. In the future this model should be extended to forest canopies by representing trees as flexible stems with various vibration modes and accounting for interactions between neighbouring trees.

We would like to thank the Center for Analysis and Prediction of Storms (CAPS) at the University of Oklahoma for providing the ARPS code. Computer simulations related to this work were performed using HPC resources from GENCI-IDRIS (Grant 2009-i2009011833) as well as Ephyse cluster. Thanks are expressed to the Ephyse computing team (Patrick Moreau, Tovo Rabemanantsoa, Guy Pracros and Dr Mark R. Irvine) for their help with the cluster set-up and administration. We would like to thank Professor Patrick Huerre for insightful discussions on mixing layers. Financial support from the ANR programme 'CHENE-ROSEAU', involving INRA, INRIA and Ecole Polytechnique, as well as support from the FQRNT, are gratefully acknowledged. Finally, we thank anonymous reviewers for their helpful comments.

## Appendix A. Subgrid-scale model

The Reynolds or subfilter-scale or subgrid stress tensor $\tau_{ij}$ is modelled through an SGS eddy-viscosity or gradient-transport model as

$$ \tau_{ij} = -((1 - \delta_{3j}) v_{th} + \delta_{3j} v_{tv}) \left( \frac{\partial \tilde{u}_i}{\partial x_j} + \frac{\partial \tilde{u}_j}{\partial x_i} \right), \quad (A 1) $$

where $v_{th}$ and $v_{tv}$ are the horizontal and vertical eddy viscosities modelled as the product of horizontal and vertical length scales, $l_{th}$ and $l_{tv}$, respectively, and a velocity scale $\sqrt{e}$ ($e$ being the SGS TKE) characterizing the SGS turbulent eddies:

$$ v_{th} = C_s \sqrt{e} l_{th}, \quad (A 2) $$

$$ v_{tv} = C_s \sqrt{e} l_{tv}, \quad (A 3) $$

where $C_s = 0.1$.
For a neutral atmosphere, the mixing lengths depend on the grid spacing:

$$l_{th} = \sqrt{\Delta x \Delta y}, \qquad (A 4)$$

$$l_{tv} = \Delta z, \qquad (A 5)$$

where $\Delta x$, $\Delta y$ and $\Delta z$ are the grid spacings in the longitudinal, lateral and vertical directions, respectively. This distinction between horizontal and vertical mixing lengths is recommended for simulations using different horizontal and vertical grid spacings.

The conservation equation for the SGS TKE writes

$$\frac{\partial e}{\partial t} + \tilde{u}_j \frac{\partial e}{\partial x_j} = -\tau_{ij} \frac{\partial \tilde{u}_i}{\partial x_j} + \frac{\partial}{\partial x_j} \left( 2 \left( (1 - \delta_{j3}) \nu_{th} + \delta_{j3} \nu_{tv} \right) \frac{\partial e}{\partial x_j} \right) \\ - \frac{g}{\theta} \tau_{3\theta} - C_{\epsilon} \frac{e^{3/2}}{l_{tv}} - 2 \frac{C_D}{l^2} \left| \tilde{u}_i - \frac{z}{h} \frac{\partial \tilde{q}_i}{\partial t} \right| e, \quad (A 6)$$

where the constant $C_{\epsilon} = 3.9$ at the lowest model level and 0.93 above accordingly to Deardorff (1980) and Moeng (1984).

The terms on the right-hand side of (A 6) represent, respectively, the dynamic shear production term, the turbulent transport term, the buoyancy production term, the dissipation term and the cascade term for SGS TKE. The latter represents the energy-loss process that accelerates the dissipation of turbulence in the canopy: as the eddies of all scales larger than the canopy elements lose their TKE into both heat and wake through their interaction with vegetation, the inertial eddy-cascade is indeed bypassed (Finnigan 2000). The production of SGS TKE by wake motions behind vegetation elements is not considered, as their scales are much smaller than those making up the bulk of SGS TKE (Shaw & Schumann 1992).

The subgrid heat flux is written as

$$\tau_{3\theta} = -\frac{\nu_{tv}}{Pr} \frac{\partial \tilde{\theta}}{\partial x_3}, \qquad (A 7)$$

where Pr is the Prandtl number.

The values of the constants $C_S$ and $C_{\epsilon}$ are entirely empirical and are derived from statistical properties of the turbulence at the SGS, which should respond to quasi isotropic eddies and a balance between the shear production and dissipation of the SGS TKE with the Kolmogorov spectrum (Moeng 1984). The values chosen for these constants in the present study are in agreement with those proposed by Deardorff (1980) and Moeng (1984). The impact of their values on the main canopy turbulent structures should be limited compared to small structures since the former structures are explicitly resolved by the model. For further details on the SGS model used by ARPS, the reader can refer to Xue et al. (2000).

## Appendix B. Discretized linear analysis equations

The z-function of the vertical velocity perturbation $\hat{w}(z)$ is discretized over the domain $[0, H]$ at $N + 2$ nodes leading to a distance between the nodes of $\Delta z = H/(N + 1)$. We can define a vector composed of the velocity coefficients, the canopy displacement and the $N$ discretized values of the z-function of the velocity perturbation away from the boundaries: $w = \langle \hat{\xi}_x, \hat{q}_x, \hat{w}_1, \hat{w}_2, \dots, \hat{w}_N \rangle^T$. We rewrite the system of equations (2.12), (2.11) and (2.10) in operator form, i.e. the eigenvalue problem of (2.13) $(A - \omega B)w = 0$. The $(N + 2) \times (N + 2)$ matrices $A$ and $B$ are the linear operators of $\tilde{w}$ corresponding to the following equations: the relation between
the coefficients of the velocity and the displacement of the canopy

$$ \hat{\xi}_x - \omega [-i\hat{q}_x] = 0, \qquad (B1) $$

the oscillator equation of the canopy motion

$$ C\hat{\xi}_x + R\hat{q}_x + \frac{i2\rho\Delta z}{kh} \sum_{n=1}^{N} \left[ \frac{d}{dz} (C_D \bar{U} z) \right]_{z_n} \hat{w}_n - i\omega M \hat{\xi}_x - \omega \frac{2\rho\Delta z}{h^2} \sum_{n=1}^{N} [C_D \bar{U} z^2]_{z_n} \hat{q}_x = 0, \qquad (B2) $$

and the governing equation of the fluid at the $j=1...N$ discretized nodes away from the boundaries

$$ \begin{aligned} & -\bar{v}_t \delta_j^4 \hat{w}_j - 2 \frac{d\bar{v}_t}{dz} \delta_j^3 \hat{w}_j + \left[ ik\bar{U} + 2k^2\bar{v}_t - \frac{d^2\bar{v}_t}{dz^2} + \frac{2}{l^2}C_D\bar{U} \right] \delta_j^2 \hat{w}_j \\ & + \left[ 4k^2 \frac{d\bar{v}_t}{dz} + \frac{2}{l^2} \frac{d}{dz} (C_D \bar{U}) \right] \delta_j \hat{w}_j + \left[ -ik^3\bar{U} - ik \frac{d^2\bar{U}}{dz^2} - k^4\bar{v}_t - k^2 \frac{d^2\bar{v}_t}{dz^2} \right] \hat{w}_j \\ & + \left[ \frac{2ik}{hl^2} \frac{d}{dz} \left( C_D \bar{U} \frac{z}{h} \right) \right] \hat{\xi}_x - \omega \delta_j^2 \hat{w}_j i + \omega \hat{w}_j ik^2 = 0, \end{aligned} \qquad (B3) $$

where the second-order-accurate centred finite differences are given by

$$ \delta_j \hat{w}_j = \frac{\hat{w}_{j+1} - \hat{w}_{j-1}}{2\Delta z}, \qquad \delta_j^2 \hat{w}_j = \frac{\hat{w}_{j+1} - 2\hat{w}_j + \hat{w}_{j-1}}{\Delta z^2}, $$

$$ \delta_j^3 \hat{w}_j = \frac{\hat{w}_{j+2} - 2\hat{w}_{j+1} + 2\hat{w}_{j-1} - \hat{w}_{j-2}}{2\Delta z^3}, $$

$$ \delta_j^4 \hat{w}_j = \frac{\hat{w}_{j+2} - 4\hat{w}_{j+1} + 6\hat{w}_j - 4\hat{w}_{j-1} - \hat{w}_{j-2}}{\Delta z^4}. $$

At the domain boundaries the conditions of (2.8) lead to

$$ \hat{w}_0 = \hat{w}_{N+1} = 0, $$

and the following ghost points which insure a second-order accuracy to the numerical scheme

$$ \hat{w}_{-1} = \frac{2}{5}\hat{w}_1 - \frac{7}{5}\hat{w}_2 + \frac{3}{5}\hat{w}_3 - \frac{1}{10}\hat{w}_4, $$

$$ \hat{w}_{N+2} = -\frac{1}{10}\hat{w}_{N-3} + \frac{3}{5}\hat{w}_{N-2} - \frac{7}{5}\hat{w}_{N-1} + \frac{2}{5}\hat{w}_N. $$

REFERENCES

ADRIAN, R. J., MEINHART, C. D. & TOMKINS, C. D. 2000 Vortex organization in the outer region of the turbulent boundary layer. *J. Fluid Mech.* **422**, 1–54.

AMIRO, B. D. 1990 Drag coefficients and turbulence spectra within 3 boreal forest canopies. *Bound.-Layer Meteorol.* **52** (3), 227–246.

ANCELIN, P., COURBAUD, B. & FOURCAUD, T. Y. 2004 Development of an individual tree-based mechanical model to predict wind damage within forest stands. *For. Ecol. Manag.* **203** (1–3), 101–121.

AUBRY, N., GUYONNET, R. & LIMA, R. 1991 Spatiotemporal analysis of complex signals: theory and applications. *J. Stat. Phys.* **64** (3–4), 683–739.

CARLOTTI, P. 2002 Two-point properties of atmospheric turbulence very close to the ground: comparison of a high resolution les with theoretical models. *Bound.-Layer Meteorol.* **104** (3), 381–410.
CASSIANI, M., KATUL, G. G. & ALBERTSON, J. D. 2008 The effects of canopy leaf area index on airflow across forest edges: large-eddy simulation and analytical results. *Bound.-Layer Meteorol.* **126** (3), 433–460.

CHOW, F. K., WEIGEL, A. P., STREET, R. L., ROTACH, M. W. & XUE, M. 2006 High resolution large-eddy simulations of flow in a steep alpine valey. Part I. Meteorology, verification and sensitivity experiments. *J. Appl. Meteorol. Climatol.* **45**, 63–86.

DEARDORFF, J. W. 1980 Stratocumulus-capped mixed layers derived from a 3-dimensional model. *Bound.-Layer Meteorol.* **18** (4), 495–527.

DOARÉ, O., MOULIA, B. & DE LANGRE, E. 2004 Effect of plant interaction on wind-induced crop motion. *Trans. ASME: J. Biomech. Engng* **126** (2), 146–151.

DROBINSKI, P., CARLOTTI, P., NEWSON, R. K., BANTA, R. M., FOSTER, R. C. & REDELSPERGER, J. L. 2004 The structure of the near-neutral atmospheric surface layer. *J. Atmos. Sci.* **61** (6), 699–714.

DROBINSKI, P., CARLOTTI, P., REDELSBERGER, J.-L., BANTA, R. M., MASSON, V. & NEWSOM, R. K. 2007 Numerical and experimental investigation of the neutral atmospheric surface layer. *J. Atmos. Sci.* **64**, 137–156.

DUPONT, S. & BRUNET, Y. 2006 Simulation of turbulent flow in an urban forested park damaged by a windstorm. *Bound.-Layer Meteorol.* **120** (1), 133–161.

DUPONT, S. & BRUNET, Y. 2008a Edge flow and canopy structure: a large-eddy simulation study. *Bound.-Layer Meteorol.* **126**, 51–71.

DUPONT, S. & BRUNET, Y. 2008b Impact of forest edge shape on tree stability: A large-eddy simulation study. *Forestry* **81**, 299–315.

DUPONT, S. & BRUNET, Y. 2008c Influence of foliar density profile on canopy flow: a large-eddy simulation study. *Agric. For. Meteorol.* **148**, 976–990.

DUPONT, S. & BRUNET, Y. 2009 Coherent structures in canopy edge flow: a large-eddy simulation study. *J. Fluid Mech.* **630**, 93–128.

DUPONT, S., BRUNET, Y. & FINNIGAN, J. J. 2008 Large-eddy simulation of turbulent flow over a forested hill: validation and coherent structure identification. *Q. J. R. Meteorol. Soc.* **134**, 1911–1929.

DWYER, M. J., PATTON, E. G. & SHAW, R. H. 1997 Turbulent kinetic energy budgets from a large-eddy simulation of airflow above and within a forest canopy. *Bound.-Layer Meteorol.* **84** (1), 23–43.

FARGE, M., SCHNEIDER, K., PELLEGRINO, G., WRAY, A. A. & ROGALLO, R. S. 2003 Coherent vortex extraction in three-dimensional homogeneous turbulence: Comparison between cvs-wavelet and pod-Fourier decompositions. *Phys. Fluids* **15** (10), 2886–2896.

FARQUHAR, T., WOOD, J. Z. & VAN BEEM, J. 2000 The kinematics of wheat struck by a wind gust. *J. Appl. Mech.* **67**, 496–502.

FESQUET, C., DUPONT, S., DROBINSKI, P., DUBOS, T. & BARTHLOTT, C. 2009 Impact of terrain heterogeneity on coherent structure properties: numerical approach. *Bound.-Layer Meteorol.* **133**, 71–92.

FINNIGAN, J. 2000 Turbulence in plant canopies. *Annu. Rev. Fluid Mech.* **32**, 519–571.

FINNIGAN, J. J. 1979 Turbulence in waving wheat. Part 1. Mean statistics and honami. *Bound.-Layer Meteorol.* **16** (2), 181–211.

FINNIGAN, J. J. & MULHEARN, P. J. 1978a Modelling waving crops in a wind tunnel. *Bound.-Layer Meteorol.* **14**, 253–277.

FINNIGAN, J. J. & MULHEARN, P. J. 1978b A simple mathematical model of airflow in waving plant canopies. *Bound.-Layer Meteorol.* **14**, 415–431.

FLESCH, T. K. & GRANT, R. H. 1991 The translation of turbulent wind energy to individual corn plant motion during senescence. *Bound.-Layer Meteorol.* **55** (1–2), 161–176.

FLESCH, T. K. & GRANT, R. H. 1992 Corn motion in the wind during senescence. Part 2. Effect of dynamic-plant characteristics. *Agronomy J.* **84** (4), 748–751.

FLESCH, T. K. & WILSON, J. D. 1999 Wind and remnant tree sway in forest cutblocks. ii. relating measured tree sway to wind statistics. *Agric. For. Meteorol.* **93** (4), 243–258.

FOSTER, R. C., VIANEY, F., DROBINSKI, P. & CARLOTTI, P. 2006 Near-surface coherent structures and the vertical momentum flux in a large-eddy simulation of the neutrally-stratified boundary layer. *Bound.-Layer Meteorol.* **120**, 229–255.
FOUDHIL, H., BRUNET, Y. & CALTAGIRONE, J. 2005 A fine-scale k-epsilon model for atmospheric flow over heterogeneous landscapes. *Environ. Fluid Mech.* **5** (3), 247–265.

GARDINER, B. A. 1995 *The Interactions of Wind and Tree Movement in Forest Canopies*. Cambridge University Press, p. 485.

GARDINER, B., PELTOLA, H. & KELLOMAKI, S. 2000 Comparison of two models for predicting the critical wind speeds required to damage coniferous trees. *Ecol. Model.* **129** (1), 1–23.

GHISALBERTI, M. & NEPF, H. 2006 The structure of the shear layer in flows over rigid and flexible canopies. *Environ. Fluid Mech.* **6** (3), 277–301.

GOSSELIN, F. & DE LANGRE, E. 2009 Destabilising effects of plant flexibility in air and aquatic vegetation canopy flows. *Eur. J. Mech. B-Fluids* **28**, 271–282.

GREEN, S. R. 1992 Modelling turbulence air flow in a stand of widely-spaced trees. *PHOENICS*, *J. Comput. Fluid Dyn.* **5**, 294–312.

HÉMON, P. & SANTI, F. 2003 Applications of biorthogonal decompositions in fluid-structure interactions. *J. Fluids Struct.* **17** (8), 1123–1143.

HOMMEMA, S. E. & ADRIAN, R. J. 2003 Packet structure of surface eddies in the atmospheric boundary layer. *Bound.-Layer Meteorol.* **106** (1), 147–170.

HUNT, J. C. R. & CARLOTTI, P. 2001 Statistical structure at the wall of the high reynolds number turbulent boundary layer. *Flow Turbul. Combust.* **66** (4), 453–475.

HUNT, J. C. R. & MORRISON, J. F. 2000 Eddy structure in turbulent boundary layers. *Eur. J. Mech. B-Fluids* **19**, 673–694.

HUTCHINS, N. & MARUSIC, I. 2007 Evidence of very long meandering features in the logarithmic region of turbulent boundary layers. *J. Fluid Mech.* **579**, 1–28.

IKEDA, S., YAMADA, T. & TODA, Y. 2001 Numerical study on turbulent flow and honami in and above flexible plant canopy. *Intl J. Heat Fluid Flow* **22** (3), 252–258.

INOUE, R. 1955 Studies of the phenomena of waving plant ('honami') caused by wind. part 1: Mechanism and characteristics of waving plant phenomena. *J. Agric. Meteorol. (Japan)* **11**, 71–82.

JAFFE, M. 1973 Thigmomorphogenesis: the response of plant growth and development to mechanical dioica. *Planta* **114**, 143–157.

JAMES, K. R., HARITOS, N. & ADES, P. K. 2006 Mechanical stability of trees under dynamic loads. *Am. J. Bot.* **93** (10), 1522–1530.

KAIMAL, J. C. & FINNIGAN, J. J. 1994 *Atmospheric Boundary Layer Flows: Their Structure and Measurements*. Oxford University Press., p. 289.

KANDA, M. & HINO, M. 1994 Organized structures in developing turbulent-flow within and above a plant canopy using a large-eddy simulation. *Bound.-Layer Meteorol.* **68** (3), 237–257.

KERZENMACHER, T. & GARDINER, B. A. 1998 A mathematical model to describe the dynamic response of a spruce tree to the wind. *Trees* **12**, 385–394.

KIM, K. C. & ADRIAN, R. J. 1999 Very large-scale motion in the outer layer. *Phys. Fluids* **11** (2), 417–422.

KRAVCHENKO, A. G., MOIN, P. & MOSER, R. 1996 Zonal embedded grids for numerical simulations of wall-bounded turbulent flows. *J. Comput. Phys.* **127** (2), 412–423.

DE LANGRE, E. 2008 Effects of wind on plant. *Annu. Rev. Fluid Mech.* **40**, 141–168.

LI, Z. & LIN, J. D. 1990 Air flow over and through a forest edge: a steady-state numerical simulation. *Bound.-Layer Meteorol.* **51**, 179–197.

LIN, C. L., MOENG, C. H., SULLIVAN, P. P. & MCWILLIAMS, J. C. 1997 The effect of surface roughness on flow structures in a neutrally stratified planetary boundary layer flow. *Phys. Fluids* **9** (11), 3235–3249.

LIU, J., CHEN, J. M., BLACK, T. A. & NOVAK, M. D. 1996 E-epsilon modelling of turbulent air flow downwind of a model forest edge. *Bound.-Layer Meteorol.* **77** (1), 21–44.

MAHRER, Y. 1984 An improved numerical approximation of the horizontal gradients in a terrain- following coordinate system. *Mon. Weather Rev.* **112**, 918–922.

MAYER, H. 1987 Wind-induced tree sways. *Trees* **1**, 195–206.

MOENG, C. H. 1984 A large-eddy-simulation model for the study of planetary boundary-layer turbulence. *J. Atmos. Sci.* **41** (13), 2052–2062.

MOENG, C.-H. & SULLIVAN, P. P. 1994 A comparison of shear-driven and buoyancy-driven planetary boundary-layer flows. *J. Atmos. Sci.* **51** (7), 999–1022.
MOULIA, B. & COMBES, D. 2004 Thigmomorphogenesis acclimation of plant to moderate winds greatly affects height structure in field-grown alfalfa (medicago sativa l.), an indeterminate herb. *Comput. Biochem. Physiol. A: Molec. Integrative Physiol.* **137** (3-S1), S77, A10.2.

PELTOLA, H. 1996 Swaying of trees in response to wind and thinning in a stand of scots pine. *Bound.-Layer Meteorol.* **77** (3-4), 285-304.

PELTOLA, H., KELLOMAKI, S., VAISANEN, H. & IKONEN, V. P. 1999 A mechanistic model for assessing the risk of wind and snow damage to single trees and stands of Scots pine, Norway spruce, and birch. *Can. J. For. Res.-Rev. Can. Rech. For.* **29** (6), 647-661.

PÉNELON, T., CALMET, I. & MIRONOV, D. V. 2001 Micrometeorological simulations over a complex terrain with submeso: a model study using a novel pre-processor. *Intl J. Environ. Pollut.* **16** (1-6), 583-602.

PY, C., DE LANGRE, E. & MOULIA, B. 2004 The mixing layer instability of wind over a flexible crop canopy. *C. R. Méc.* **332** (8), 613-618.

PY, C., DE LANGRE, E. & MOULIA, B. 2006 A frequency lock-in mechanism in the interaction between wind and crop canopies. *J. Fluid Mech.* **568**, 425-449.

PY, C., DE LANGRE, E., MOULIA, B. & HEMON, P. 2005 Measurement of wind-induced motion of crop canopies from digital video images. *Agric. For. Meteorol.* **130** (3-4), 223-236.

RAUPACH, M. R., FINNIGAN, J. J. & BRUNET, Y. 1996 Coherent eddies and turbulence in vegetation canopies: the mixing-layer analogy. *Bound.-Layer Meteorol.* **78** (3-4), 351-382.

RODRIGUEZ, M., DE LANGRE, E. & MOULIA, B. 2008 A scaling law for the effects of architecture and allometry on tree vibration modes suggests a biological tuning to modal compartmentalization. *Am. J. Bot.* **95** (12), 1523-1537.

ROSS, A. N. 2008 Large-eddy simulations of flow over forested ridges. *Bound.-Layer Meteorol.* **128**, 59-76.

RUSSELL, G., MARSHALL, B. & JARVIS, P. G. 1990 *Plant Canopies: Their Growth, Form and Function*. Cambridge University Press, p. 192.

SELLIER, D., BRUNET, Y. & FOURCAUD, T. 2008 A numerical model of tree aerodynamic response to a turbulent airflow. *Forestry* **81** (3), 279-297.

SELLIER, D., FOURCAUD, T. & LAC, P. 2006 A finite element model for investigating effects of aerial architecture on tree oscillations. *Tree Physiol.* **26** (6), 799-806.

SHAW, R. H., BRUNET, Y., FINNIGAN, J. J. & RAUPACH, M. R. 1995 A wind tunnel study of air flow in waving wheat: two-point velocity statistics. *Bound.-Layer Meteorol.* **76** (4), 349-376.

SHAW, R. H. & SCHUMANN, U. 1992 Large-eddy simulation of turbulent-flow above and within a forest. *Bound.-Layer Meteorol.* **61** (1-2), 47-64.

SHEN, S. H. & LECLERC, M. Y. 1997 Modelling the turbulence structure in the canopy layer. *Agric. For. Meteorol.* **87** (1), 3-25.

SU, H. B., SHAW, R. H. & PAW, U, K. T. 2000 Two-point correlation analysis of neutrally stratified flow within and above a forest from large-eddy simulation. *Bound.-Layer Meteorol.* **94** (3), 423-460.

SU, H.-B., SHAW, R. H., PAW, U, K. T., MOENG, C.-H. & SULLIVAN, P. P. 1998 Turbulent statistics of neutrally stratified flow within and above a sparse forest from large-eddy simulation and field observations. *Bound.-Layer Meteorol.* **88**, 363-397.

WATANABE, T. 2004 Large-eddy simulation of coherent turbulence structures associated with scalar ramps over plant canopies. *Bound.-Layer Meteorol.* **112** (2), 307-341.

WEIGEL, A. P., CHOW, F. K. & ROTACH, M. W. 2007 On the nature of turbulent kinetic energy in a steep and narrow alpine valley. *Bound.-Layer Meteorol.* **123** (1), 177-199.

WEIGEL, A. P., CHOW, F. K., ROTACH, M. W., STREET, R. L. & XUE, M. 2006 High-resolution large-eddy simulations of flow in a steep alpine valley. Part II. Flow structure and heat budgets. *J. Appl. Meteorol. Climatol.* **45** (1), 87-107.

XUE, M., DROEGEMEIER, K. K. & WONG, V. 2000 The advanced regional prediction system (arps) – a multi-scale nonhydrostatic atmospheric simulation and prediction model. Part I. Model dynamics and verification. *Meteorol. Atmos. Phys.* **75** (3-4), 161-193.

XUE, M., DROEGEMEIER, K. K., WONG, V., SHAPIRO, A. & BREWSTER, K. 1995 *ARPS version 4.0 Users Guide*. Center for Analysis and Prediction of Storms, p. 380.

XUE, M., DROEGEMEIER, K. K., WONG, V., SHAPIRO, A., BREWSTER, K., CARR, F., WEBER, D., LIU, Y. & WANG, D. 2001 The advanced regional prediction system (arps) – a multi-scale nonhydrostatic
atmospheric simulation and prediction tool. Part II. Model physics and applications. *Meteorol. Atmos. Phys.* **76** (3–4), 143–165.

YANG, B., MORSE, A. P., SHAW, R. H. & PAW, U. K. T. 2006a Large-eddy simulation of turbulent flow across a forest edge. Part II. Momentum and turbulent kinetic energy budgets. *Bound.-Layer Meteorol.* **121** (3), 433–457.

YANG, B., RAUPACH, M. R., SHAW, R. H., THA, K., PAW, U. & MORSE, A. P. 2006b Large-eddy simulation of turbulent flow across a forest edge. Part I. Flow statistics. *Bound.-Layer Meteorol.* **120** (3), 377–412.