Boundary-Layer Meteorology (2020) 176:35–60
https://doi.org/10.1007/s10546-020-00519-8
RESEARCH ARTICLE

Large-Eddy Simulation of Neutrally-Stratified Turbulent Flow
Within and Above Plant Canopy Using
the Central-Moments-Based Lattice Boltzmann Method
Tsutomu Watanabe1 · Kou Shimoyama1 · Masayuki Kawashima1 ·Yasuko Mizoguchi2 · Atsushi Inagaki3
Received: 27 September 2019 / Accepted: 24 March 2020 / Published online: 13 April 2020
© The Author(s) 2020

Abstract
A central-moments-based lattice Boltzmann model for large-eddy simulation of neutrallystratified turbulent flows is described. Through comparative simulations of the airflow within
and above a homogeneous plant canopy, the performance of the model is evaluated withrespect to a conventional large-eddy-simulation model based on the incompressible Navier–Stokes equations. Simulated turbulence statistics, such as the mean velocity, velocity variances,velocity skewness, and power spectra, are shown to be almost identical between the twomodels. The spatial structure of coherent eddies and their maintenance processes are also
confirmed to be properly represented by the lattice Boltzmann method through analysis of
the turbulence kinetic energy budget and spatial two-point correlation functions. Using thesimulated results, the energetics of the streamwise-elongated streaky structures commonly
observed over vegetation and urban canopies are examined. While the short-wavelengthcomponents of the shear-generated streamwise kinetic energy are redirected rapidly by pressure to the lateral and vertical velocity components, long-wavelength energy tends to remainin the streamwise velocity component, which is dissipated in relatively slower processes.
Consequently, the streaky structures persist in the streamwise velocity component.
Keywords Coherent structure · Large-eddy simulation · Lattice Boltzmann method · Plant
canopy · Turbulence

B Tsutomu Watanabe

t-wata@lowtem.hokudai.ac.jp1

Institute of Low Temperature Science, Hokkaido University, Kita 19 Nishi 8, Kita-ku, Sapporo
060-0819, Japan

2

Hokkaido Research Center, Forestry and Forest Products Research Institute, 7 Hitsujigaoka,
Toyohira, Sapporo 062-8516, Japan3

Department of Transdisciplinary Science and Engineering, Tokyo Institute of Technology, 2-12-1
O-okayama, Meguro-ku, Tokyo 152-8552, Japan

123

36

T. Watanabe et al.

1 Introduction
The exchange of momentum, heat, water vapour, carbon dioxide and other gasses, and particles between vegetation and the atmosphere, is one of the main processes in land–atmosphereinteractions. Turbulent eddies promote vegetation–atmosphere exchange by transferring thesequantities through the plant canopy layer, in which individual sources and sinks are spatiallydistributed below plant height, to and from the roughness sublayer above the canopy.
In neutrally-stratified conditions, turbulence both in the roughness sublayer and in the
canopy layer originally stems from strong mean shear, which is the consequence of the
reduction of wind speed in the canopy layer due to the drag force exerted by canopy elements. Turbulence kinetic energy (TKE) produced near the canopy top is then transported by
eddies into the canopy layer, where it is dissipated substantially in the wake of canopy elements (Raupach and Thom 1981). The spatial length scale of most energetic shear-generated
eddies is typically a few times larger than the canopy height (Raupach et al. 1996). Thesephenomenological characteristics explain a number of well-documented features of canopy
turbulence that differ from those of ordinary boundary-layer turbulence. Such features include
large values of skewness in the streamwise and vertical velocity components, which are consistent with the dominance of sweep motions near the canopy top, velocity flatness factorsthat greatly exceed the Gaussian value, the marked contribution of vertical transport in thebudgets of second-order moments, and counter-gradient flux phenomena, indicating the inadequacy of the local gradient-diffusion model for representing the fluxes of momentum andscalars in the canopy layer (Shaw 1977; Wilson and Shaw 1977; Raupach and Thom 1981;
Denmead and Bradley 1985; Watanabe 1993; Raupach et al. 1996; Katul et al. 2013). All
these features are consistent with the fact that tall, coherent eddies make essential contributions to the vertical fluxes (Denmead and Bradley 1985; Bergstöm and Högström 1989; Gaoet al. 1989), and constitute the majority of turbulent structures within and above the canopy
layer (Shaw et al. 1995; Finnigan and Shaw 2000).
In addition to field observations, flume or wind-tunnel experiments, and theoretical considerations, numerical analysis using large-eddy simulation (LES) has received increasing
attention in the previous three decades. Since Shaw and Schumann (1992) and Kanda and
Hino (1994) first conducted LES investigations of turbulent flows within and above plant
canopies, many other authors have published reports on numerical analyses of canopy turbulence in homogeneous plant canopies (e.g., Dwyer et al. 1997; Shen and Leclerc 1997; Suet al. 1998; Patton et al. 1998; Fitzmaurice et al. 2004; Watanabe 2004, 2009; Yue et al. 2007;
Finnigan et al. 2009; Shaw et al. 2013), heterogeneous plant canopies (Dupont and Brunet
2009; Bohrer et al. 2009), and a homogeneous plant canopy under the influence of convective
motions in the atmospheric boundary layer (Patton et al. 2016). The LES approach has also
been applied to urban canopies in pioneering work by Kanda et al. (2004), Kanda (2006), andCoceal et al. (2006), followed by Inagaki et al. (2012) and others, and LES is now established
as an indispensable tool for investigation of the nature of canopy turbulence.
The LES models used in previous studies have generally been based on the Navier–Stokes
equations and the incompressible continuity equation, which is indeed the most orthodox
method in computational fluid dynamics. As these equations directly express macroscopic
conservation laws, it is straightforward to compile simulation code that conserves momentum
and mass. Kinetic energy can also be conserved if the non-linear advection terms are treatedproperly. Therefore, such models are superior in terms of property conservation, which isdesirable for the reliable simulation of turbulence.

123Large-Eddy Simulation of Neutrally-Stratified Turbulent Flow37

The lattice Boltzmann method (hereafter, the "Boltzmann approach") is a recently
advanced method of computational fluid dynamics (e.g., Chen and Doolen 1998; Aidum andClausen 2010), which, instead of considering macroscopic velocities and density, describes
fluid motions by the spatiotemporal variations of the discrete velocity-distribution function of
microscopic fluid particles (Krüger et al. 2017). In an Eulerian framework, the temporal evolution of the distribution function at a given node in a regular grid (lattice) system is governedby the discretized Boltzmann equation, which incorporates particle streaming, collision, andforcing processes as the key physics. Macroscopic variables (velocities and density) are evaluated as low-order moments of the distribution function. It can be shown mathematically thatthe derived macroscopic variables satisfy the weakly compressible Navier–Stokes equations
within the limit of a low Mach number (Chen and Doolen 1998; Krüger et al. 2017). The
lattice Boltzmann method has several advantages over the incompressible Navier–Stokes
models. Solving the Poisson equation is not necessary for evaluation of the pressure, whichis calculated simply from an equation of state, i.e., pressure is related linearly to density
in an isothermal condition. Advection of macroscopic variables can be simulated with highaccuracy because the streaming process in the method is expressed exactly as a translationof the distribution function during one timestep to a neighbouring grid point in the direction of the prescribed constant microscopic velocities. Collision and forcing processes are
expressed in a local manner, in which reference is only to a set of variables defined at asingle grid point. With these features, together with its simple algorithm and the ease of
implementing boundary conditions, the lattice Boltzmann method is feasible for application
to high-resolution simulations of atmospheric-boundary-layer flows using advanced parallelprocessing machines.
Application of the lattice Boltzmann method to the atmospheric-boundary-layer flow was
initiated only recently. Onodera et al. (2013) compiled a LES code by incorporating an
eddy-viscosity model (Kobayashi 2005) into a lattice Boltzmann algorithm, and conducteda precise numerical simulation of airflow within and above an actual urban canopy. Throughusing a massively parallelized supercomputer, the authors simulated the detailed wind field
at 1-m resolution over individual buildings in an area of approximately 10 km × 10 km in
the city of Tokyo, Japan. Ahmad et al. (2017) and Inagaki et al. (2017) used the same code
to perform simulations of the spatial development of a neutral boundary layer over an actual
urban geometry near Tokyo with 2-m resolution in a domain of 19.2 × 4.8 × 1.0 km3 . The
latter simulation allowed the authors to conduct simultaneous analysis of the aerodynamic
properties both of boundary-layer-scale eddies and of gusty motions in street canyons. Thesestudies demonstrated the potential of the lattice Boltzmann approach to be a powerful tool
for simulations of turbulent flows in the atmospheric boundary layer.
This paper describes the development of a new LES code based on the lattice Boltzmann
method for the simulation of canopy turbulence in neutrally-stratified conditions, presents
the validation of the code with respect to a conventional LES model based on the Navier–Stokes equations, and briefly discusses the maintenance process of the large-scale eddies that
dominate canopy turbulence.

2 Model Description
We adopt a Cartesian coordinate system in which x denotes the streamwise direction, yrepresents the lateral direction, and z represents the vertical direction. Velocity components

123

38

T. Watanabe et al.

in these three directions are u, v, and w, respectively. Tensor notation such as (x 1 , x2 , x3 ) (x, y, z) and (u 1 , u 2 , u 3 )  (u, v, w) is also used.

2.1 Navier–Stokes-Based Model
A LES model described by Watanabe (2004, 2009) is used for the reference simulation, andis based on the Navier–Stokes equations and the incompressible continuity equation. The
subgrid-scale (SGS) TKE is modelled using a separate prognostic equation. The SGS TKE
is used to model the eddy viscosity, which relates the strain-rate tensor of resolved velocities
to the SGS momentum fluxes. Horizontal derivatives are represented by a pseudo-spectral
method with the 2/3 dealiasing rule, while a fourth-order Padé scheme is adopted for verticaldifferencing. Time advancement is performed by the third-order Runge–Kutta scheme. At
each step in the Runge–Kutta scheme, the Poisson equation for pressure is solved in Fourier
space to satisfy the constraint of incompressible continuity. The plant canopy is modelled
as homogeneously distributed sinks for momentum and SGS TKE. The instantaneous sinkstrengths are described by the plant area density, an isotropic drag coefficient, and the resolved
velocity components. The numerical set-up for the reference simulation is the same as thatof the lattice Boltzmann method (described later), except that the exact value of the external
driving force is different and a longer streamwise domain (1500, 400, and 81 nodes in thestreamwise, lateral, and vertical directions, respectively) is adopted because this simulationis conducted for the further purpose of inspecting the longest streamwise modes of thevelocity fluctuations. However, an appropriate normalization described in Sect. 3 enables aproper comparison of turbulence statistics presented here. The influence of the domain size
is reported elsewhere (Watanabe 2004, 2009).

2.2 Lattice Boltzmann Method
General and detailed descriptions of the lattice Boltzmann method can be found in publishedreviews (e.g., Chen and Doolen 1998; Aidum and Clausen 2010), textbooks (e.g., Mohamad
2011; Krüger et al. 2017), and references therein. Only the essence adequate for reproducing
the present model is given here.
The basic structure of our lattice Boltzmann method is a three-dimensional 27-velocity(D3Q27) lattice model, which is considered suitable for simulations of high Reynolds number flows in terms of numerical stability and accuracy (Suga et al. 2015). In this model,
the distribution functions f i jk (x, y, z, t) at position x and time t of discrete microscopic
velocities (ic, jc, kc) are considered. The indices i, j, k ∈ {−1, 0, 1} represent the direction of the discretized velocities (thus 27 velocities in total), and c  x/t is the lattice
velocity unit defined with lattice spacing x and timestep t. The main part of the lattice
Boltzmann algorithm consists of three steps: collision, forcing, and streaming. The collision
step represents the tendency of distribution functions towards equilibrium (with respect tomacroscopic velocities and density) due to collisions between particles. The forcing step calculates temporal variations of distribution functions induced by external and internal forces.
The streaming step describes the translation of fluid particles in the direction of the discrete
microscopic velocities. The equations involved in each step are shown below. Following theconvention in the literature, the equations are written using the lattice-unit system in which
all variables are normalized, such that x  t  c  1 should result and the mean fluid
density should be equal to unity. Therefore, the differences x and t are eliminated frommost of the following equations for brevity.

123

Large-Eddy Simulation of Neutrally-Stratified Turbulent Flow39

The collision process is modelled using a central-moment-based multi-relaxation-timemodel (Geier et al. 2006, 2015; Premnath and Banerjee 2009, 2011), which is known tohave higher numerical accuracy and stability compared with more commonly used singlerelaxation-time models and other multi-relaxation-time models based on raw moments (Geieret al. 2015). A cumulant-based multi-relaxation-time model (Geier et al. 2015; 2017) isanother possible candidate if a few more computations can be devoted for the collision step.
The central moments of the distribution function, κ pqr , are defined asκ pqr (x, t) 

1
1
1




(ic − u) p ( jc − v)q (kc − w)r f i jk (x, t),(1)i−1 j−1 k−1where the indices p, q, r ∈ {0, 1, 2} are the directional orders of moment, and the sum of theindices ( p + q + r ) defines the formal order of moment (e.g., κ200 , κ110 , and other momentswith permutated indices are of the second order).
In the lattice Boltzmann algorithm, the distribution functions and the macroscopic variables (density, velocities, and pressure) are updated through the following procedure. At
every timestep, the distribution functions are transformed by Eq. 1 into the central-moment
space, where the collision and forcing steps proceed asκ ∗pqr (x, t)  κ pqr (x, t) + κpqr (x, t) + S κpqr (x, t).

(2)Here, the asterisk indicates the state after the collision and forcing steps (hereafter, the“post-collision state”), and κpqr and S κpqr respectively denote the collision and forcing
terms described later. The post-collision central moments are then transformed back to the
distribution function space byf ∗ (x, t)  C −1 κ∗ (x, t),(3)where C −1 represents the inverse transform of Eq. 1. In the streaming step, the post-collision
distribution functions are translated to adjacent grid points as
f i jk (x + ict, y + jct, z + kct, t + t)  f i∗jk (x, y, z, t).

(4)

After the streaming step, which is followed by application of appropriate boundary conditions,the macroscopic variables at the next timestep are calculated from the low-order raw momentsof the distribution function as
ρ
ρu 
ρv 
ρw 

1
1
1




f i jk ,(5)(ic) f i jk + Fx /2,(6a)( jc) f i jk + Fy /2,(6b)(kc) f i jk + Fz /2,(6c)i−1 j−1 k−1
1
1
1



i−1 j−1 k−1
1
1
1



i−1 j−1 k−1
1
1
1



i−1 j−1 k−1p  ρ/3,(7)

where p is the macroscopic static pressure.

123

40

T. Watanabe et al.

As explained, the collision step describes the tendency of the velocity distribution towards
equilibrium. The continuous Maxwell–Boltzmann distribution f M defined as 3/2 
3
3|ξ − u|2M
exp −f (ρ, u, ξ)  ρ
,
(8)
2π
2


where ξ  ξx , ξ y , ξz denotes the microscopic particle velocity, is adopted as the equilibrium
distribution to which the velocity distribution is relaxed. The continuous central moments ofeq
the equilibrium distribution, κ pqr , are defined aseqκ pqr +∞˚q(ξx − u) p ξ y − v (ξz − w)r f M dξx dξ y dξz ,

(9)

−∞

and are shown analytically to have non-zero values only for the moments for which thedirectional orders (p, q, r) are all even (Geier et al. 2006; Premnath and Banerjee 2011), i.e.,eq

κ000  ρ,
eq

eq

eq

eq

eq

eq

κ200  κ020  κ002  ρ/3,
κ220  κ022  κ202  ρ/9,eqκ222  ρ/27.

(10)

The forcing processes are also considered in the central-moment space. Defining F Fx , Fy , Fz as a force vector, the temporal rate of the distribution function attributable
to that force can be approximated as f /t  (3/ρ) (ξx − u)Fx + ξ y − v Fy + (ξz − w)Fz f M ,(11)based on the assumption that the particle-velocity distribution is generally near the equilibrium state represented by the Maxwell–Boltzmann distribution (He et al. 1998b; Premnathand Banerjee 2009). The forcing term S κpqr is given the continuous central moments of Eq. 11,
which can also be derived analytically and the non-zero values are written asκκκ Fx , S010
 Fy , S001
 Fz ,
S100
κ
κ
S120
 S102
 Fx /3,

κκ
S210
 S012
 Fy /3,
κκ
S021
 S201
 Fz /3,

κ
κ
κ
S122
 Fx /9, S212
 Fy /9, S221
 Fz /9.

Using these expressions, the collision term is basically modelled as


1 κ
eq
κ
 pqr  −ω pqr κ pqr + S pqr − κ pqr ,
2

(12)

(13)

where ω pqr denotes the collisional relaxation coefficients for the individual moments. Equation 13 represents that the central moments after the half forcing is relaxed towards theeqequilibrium moments κ pqr , which are associated with the macroscopic density and velocities(Eqs. 8 and 9) that are also affected by the same amount of forcing, as implied by Eq. 6. By
adopting this representation, the overall numerical accuracy can be maintained at the second
order in both space and time (Guo et al. 2002; Krüger et al. 2017). However, as applying Eq. 13
independently to each of the moments may break the rotational invariance, some moments
are linearly combined and classified into 10 groups for which the same relaxation coefficients123

Large-Eddy Simulation of Neutrally-Stratified Turbulent Flow41

(ω1 , ω2 , . . . , ω10 ) are assigned (Geier et al. 2015). The equations actually implemented forthe calculation of the post-collision moments are fully described in Appendix.

2.3 Subgrid-Scale Viscosity and Relaxation Coefficients
The macroscopic velocities predicted by the lattice Boltzmann method are defined at discretenodes of the lattice system and therefore correspond to the resolved-scale components. The
SGS stress is modelled by the eddy-viscosity model as∂u α ∂u β2
ταβ  −υSGS+ kSGS δαβ ,+(14)∂ xβ ∂ xα3where α and β are indices for spatial coordinates, ταβ is the SGS stress, υSGS is the SGS eddy
viscosity, kSGS is the SGS kinetic energy, and δαβ is the Kronecker delta. A Smagorinsky-type
parametrization is used for the subgrid eddy viscosity,
υSGS  CSGS 2 2Sαβ Sαβ ,(15)where CSGS is a parameter that depends on time and space,  is the grid spacing, andSαβ is the resolved strain-rate tensor; the repetition of indices implies summation over all
three components unless otherwise stated. The strain-rate tensor in Eq. 15 is evaluated using
local information (Eqs. 48 and 49 in Appendix) rather than finite differences that refer to themacroscopic velocities at neighbouring grid points. In accordance with the coherent-structure
Smagorinsky model proposed by Kobayashi (2005), the parameter CSGS is evaluated by
CSGS  C1

Q 3/2
,
E

(16)

√
where C 1  3/20 is a model parameter optimized to reproduce the velocity spectra well
above the canopy in preliminary simulations, Q is the second invariant of the velocity gradienttensor, and E is the magnitude of the velocity gradient tensor,
1
Q  Wαβ Wαβ − Sαβ Sαβ ,
(17)
2

1
(18)
E  Wαβ Wαβ + Sαβ Sαβ ,2with Wαβ being the resolved vorticity tensor. Only for the evaluation of the tensors Q and
E in these equations are the velocity gradients estimated using the second-order centred
difference. As the value of CSGS given by Eq. 16 automatically vanishes on a solid wall
(Kobayashi 2005), wall-damping functions are not introduced.
In the framework of the lattice Boltzmann approach, the SGS eddy viscosity is incorporated
into the relaxation coefficient ω1 for a group of the second-order moments (see Appendix)
as (Aidum and Clausen 2010)


1 1
1
υ + υSGS 
,
(19)
−
3 ω1
2
where υ is the kinematic shear viscosity. Relaxation coefficients for other groups of the
moments are practically free parameters in the range of 0 < ω < 2, because they do notappear in the Navier–Stokes equations recovered from the lattice Boltzmann equations (seebelow) and thus have no primary influence on the simulated velocity fields. The values usedhere are ω2  ω6  ω7  ω8  ω10  ω1 and ω3  ω4  ω5  ω9  1, which
are determined through preliminary simulations to achieve suitable numerical stability as

123

42

T. Watanabe et al.

simple as possible. It should be noted that, although the parameter υ is included in Eq. 19
for generality, the effect of shear viscosity is negligible in the present simulation of a welldeveloped turbulent flow with a Reynolds number of Re ~ 105 (based on the canopy height
and the friction velocity defined later).
To account correctly for the second term on the right-hand side of Eq. 14, the SGS kinetic
energy kSGS must be evaluated. Following Inagaki et al. (2010) and Suga et al. (2015), thevalue of kSGS is estimated by
kSGS  ρck

3
(u α − u α )2 ,

(20)

α1

where ck ( 1) is a model constant and u α  is a filtered velocity component given by

1 B
E
N
T
u + u Sα + u W
u α  (21)
α + 6u α + u α + u α + u α ,12 αSin which u B
α , u α , . . . denote velocities at the six surrounding grid points. By adding the term−(2/3)∂kSGS /∂ xα to all the force components (Fα ) in the above equations, the effects of thenormal components of the SGS stress are included properly.
With these parametrizations, the present lattice Boltzmann model predicts the macroscopic
density and velocities that satisfy the filtered continuity equation,∂ρ ∂ρu α
+
 0,
∂t
∂ xα

(22)

and the filtered Navier–Stokes equations,∂ρu α ∂ρu α u β∂u α ∂u β2∂p∂ρ(υ + υSGS )− kSGS δαβ + Fα . (23)
+
−
+
+
∂t
∂ xβ
∂ xα ∂ xβ∂ xβ ∂ xα
3

2.4 Set-Up for the Canopy-Turbulence Simulation
The lattice Boltzmann simulation is performed for a turbulent flow within and above ahomogeneous plant canopy in a computational domain of 1024, 512, and 80 nodes in thestreamwise, lateral, and vertical directions, respectively. The lowest 10 nodes comprise thecanopy layer, in which the effect of plant elements is considered a spatially distributed
momentum sink. The grid spacing therefore corresponds to 0.1 h, with h being the canopy
height.
The airflow in the computational
domain is driven by a streamwise external force F ext ,
extwhich is fixed constant at F / ρx/(t)2  10−7 . Horizontal boundaries are periodic,
while the top boundary is a rigid free-slip wall. The bottom boundary (i.e., the ground surface)
is a rigid wall with friction that is represented by including an instantaneous frictional force
in the horizontal force terms at the lowest fluid nodes. The rigid (and free-slip) boundaries are
represented by means of halfway specular reflection, in which particle-velocity componentsthat move inwards from the wall are given by
f i jk (x, y, z b , t + t)  f i∗jk (x − ic, y − jc, z b , t),(24)

where z b is the height of the highest or lowest fluid node, and k denotes an inwards direction
opposite to the direction towards the wall (k ). This implies that the actual wall boundaries
are located one-half grid spacing beyond the highest or lowest fluid node (Lallemand andLuo 2003); all variables ( f , u, v, w,ρ, p) are collocated at each grid node.

123

Large-Eddy Simulation of Neutrally-Stratified Turbulent Flow43

The instantaneous drag force exerted by the plant elements is conventionally modelledusing a product of the leaf area density a, isotropic drag coefficient cd , and instantaneous
macroscopic velocity u as
Fd  −ρcd a|u|u.

(25)

These canopy parameters are given as ah  2 (i.e., the leaf area index LAI  2) and cd  0.2.
Similarly, the frictional force of the underlying ground is given by the bulk relation,
Fg  −ρCM |ub |ub ,(26)

where CM is the bulk transfer coefficient of momentum and ub is the instantaneous macroscopic velocity at the lowest nodes. The coefficient CM is calculated assuming a logarithmic
wind profile with a prescribed aerodynamic roughness length z 0g , which is prescribed as
z 0g /h  10−3 throughout the bottom boundary. The drag force depends on the square of themacroscopic velocity (Eq. 25). However, as the definition of macroscopic velocity includes
the drag force (Eq. 6), the macroscopic velocity is evaluated as a solution of the resulting
quadratic equation (Guo and Zhao 2002),
u0

u
,
(27)
1
1
1
+
+
c
a|u
|
d
0
2
4
2
where u0 denotes the macroscopic velocity calculated without the drag force from Eq. 6. At
the lowest nodes, cd a in the above equation is replaced by (cd a + CM ) to account for thefrictional force of the ground.
Each simulation is initialized with a uniform streamwise velocity component with a small
amount of random perturbations in the vertical velocity component and the distributionfunctions that represent equilibrium to the given velocities. After the turbulence statistics
reach a steady state, three-dimensional data of the macroscopic variables and the distributionfunctions are collected at a constant time interval of 0.66 h/u ∗ . In total, 450 samples are used
for the analysis described in the next section.

2.5 Implementation of the Lattice Boltzmann Method
In the actual implementation, the distribution functions pre-conditioned by subtracting constant values are used as suggested by Geier et al. (2015, 2017) to minimize round-off errors.
Moreover, the corresponding transformations to/from the central-moment space are performed using a factorized fast procedure described by Geier et al. (2017). The model, codedusing the Compute Unified Device Architecture (CUDA) C language, runs on a desktop
PC (CPU: Intel Core i7-4790K) equipped with a graphics processing unit (NVIDIA Quadro
GV100). Double-precision (64 bit) arithmetic is used. As suggested by Rinaldi et al. (2012)and Delbosc et al. (2014), a streaming–collision combined kernel is adopted to reduce thenumber of accesses to the global memory and thereby to improve overall computational
efficiency.

3 Results and Discussions
The results are presented non-dimensionalized using the friction velocity u ∗ and canopy
height h. In a horizontally homogeneous steady airflow driven by a constant volumetricforce, the mean shear stress above the canopy should decrease linearly with height to zero at

123

44

T. Watanabe et al.

the free-slip top boundary. The shear stress at the canopy top τh can be predicted theoreticallyasτh  F ext z top − h ,(28)where z top is the top-boundary√height ( 8h). The friction velocity in the present simulations
is therefore defined by u ∗  τh /ρ. In what follows, an overbar ( ) denotes a mean value
obtained by averaging over horizontal planes at a given height in all samples, while a prime( ) indicates deviations from the mean.

3.1 One-Dimensional Profile
Vertical profiles of the mean wind speed, turbulent momentum flux, resolved velocity variances, and resolved velocity skewness are shown in Fig. 1. The mean wind speed decreases
exponentially with depth into the canopy and increases near-logarithmically with height
above the canopy. An inflection point of the profile is thus located near the canopy top.
As expected, the magnitude of the vertical momentum flux decreases linearly with height
above the canopy, while it decreases quickly within the canopy because of the existenceof a momentum sink. Velocity variances have maxima at z/h  1.2–1.8 and their relative
magnitudes are u u > v v > w w . Velocity skewness in the canopy layer is positive forthe streamwise velocity component and negative for the vertical velocity component, andthe signs become reversed at a distance well above the canopy. The skewness of the lateralvelocity component is essentially zero at all heights. All these characteristics share common
well-known features of turbulent flows within and above plant canopies (Kaimal and Finnigan 1994; Raupach et al. 1996). While the profiles for these quantities are similar to thosefrom the Navier–Stokes model, a remarkable difference can be seen in Fig. 1c, in which
the Navier–Stokes u u profile exhibits an oscillation near the canopy top. Because of the
stepwise change given in the leaf area density, the velocity varies sharply across the canopy
top. In such a circumstance, vertical advection near the canopy top is approximated lessaccurately by the finite-difference method adopted in this Navier–Stokes model (Watanabe
2004, 2009) and dispersive errors grow around this point; a similar result was found by Ouwersloot et al. (2017). To reduce such errors, a higher-order less dispersive scheme must be
used. Conversely, streaming of the distribution functions in the lattice Boltzmann method
is always performed by an exact upstream scheme with a Courant number of unity (Eq. 4),
which enables the model to simulate three-dimensional advection with almost no numericalerror. In general, this is one of the remarkable advantages of the lattice Boltzmann methodin simulating airflow within and above canopies with complex geometries.

3.2 Turbulence Kinetic Energy Budget
Figure 2 shows vertical profiles of the individual terms in the budget equation of the TKE,∂u∂ρw e ∂ρu α τ α3
∂p w
∂ρe
0−
−
−
− ρu w− cd aρ|u|u α u α − ρε,
∂t
∂z
∂z
∂z
∂z

(29)

where ρe denotes the resolved TKE with e  u α u α /2, and ρε is the loss of TKE due to the
production of SGS motion. The terms on the right-hand side of Eq. 29 represent (from left toright) turbulent transport, SGS transport, pressure transport, shear production, loss of TKE
due to the production of small-scale turbulence in the wakes of canopy elements, and SGS
production, respectively. The vertical fluxes in this equation are calculated at the same grid123Large-Eddy Simulation of Neutrally-Stratified Turbulent Flow

45

Fig. 1 Vertical profiles of a mean wind speed, b momentum flux, c resolved velocity variances, and d resolved
velocity skewness calculated from the lattice Boltzmann method (LBM) (mark and line) and the Navier–Stokes(NS) model (line). Broken lines in the (b) represent the SGS components of the momentum flux, and grey
shading indicates the canopy layer in which the leaf area is distributed homogeneously

points as the velocities, while the vertical gradients are evaluated at the grid boundaries. The
SGS production term for the Navier–Stokes model is calculated directly fromε  τ αβ ∂u α /∂ xβ ,(30)

whereas the same term for the lattice Boltzmann method is evaluated as the residual of
Eq. 29. In Fig. 2, the Navier–Stokes results are smoothed by a vertical 1–2–1 filter to reduce

123

46

T. Watanabe et al.

Fig. 2 Vertical profiles of the terms in the TKE budget equation simulated by the lattice Boltzmann (lines) andthe Navier–Stokes (marks) models. Ps : shear production, Tt : turbulent transport, Tp : pressure transport, Dc :loss by canopy drag, TSGS : SGS transport, DSGS : SGS production. A non-zero Dc value for the Navier–Stokes
model above the canopy is attributable to the 1–2–1 filtering

numerical oscillations and to facilitate comparison with the lattice Boltzmann results. The
figure implies that TKE is injected by the mean shear and consumed by the SGS productions
above the canopy and the wake production within the canopy. Vertical transfer by the resolved
turbulence, SGS turbulence, and pressure work in tandem to redistribute TKE to deep inside
the canopy layer. Although these features are well known in relation to canopy turbulence(Wilson and Shaw 1977; Raupach et al. 1986; Meyers and Baldocchi 1991; Dwyer et al.
1997), this figure clearly indicates that the method is as capable as the Navier–Stokes modelof reproducing the key processes of production, dissipation, and transport of TKE.
Certain differences between the two models can be seen in terms of the pressure transportand SGS production near the canopy top, with the difference in pressure transport presumably
attributable to the differences of the model equations. The Navier–Stokes model is based on
the incompressible continuity equation, whereas the lattice Boltzmann equation possesses a
weakly compressible nature (Eqs. 22 and 23), which permits the propagation of pseudo soundwaves. The reason why the difference in SGS production cannot be determined explicitly isthat this term in the lattice Boltzmann equation is a residual of the TKE budget. One possible
factor might be the aforementioned dispersion error associated with the Navier–Stokes model
used here. If so, the numerical dispersion locally increases the spatial derivatives of velocities
near the canopy top, thereby producing a local peak of SGS production at that point.

3.3 Spectra
Streamwise and lateral one-dimensional power spectra of the velocity components u, v, andw at selected levels are shown in Fig. 3. Spectra at levels above the canopy (z/h  3.05,
2.05) generally follow the −5/3 power law in the inertial subrange, whereas spectra within
the canopy (z/h  0.55) have steeper gradients due to the so-called short-circuiting processof the inertial energy cascade caused by the drag of canopy elements (Raupach and Shaw
1982; Finnigan 2000). As obstacles such as stems and branches of a particular length scale123

Large-Eddy Simulation of Neutrally-Stratified Turbulent Flow47

are not represented in the model, a secondary peak in the inertial subrange that is sometimes
observed inside canopy layers (Seginer et al. 1976; Raupach et al. 1986; Poggi et al 2004;
Cava and Katul 2008; Dupont et al. 2012) is not seen in Fig. 3. High-performance simulations
that resolve small-scale eddies shed by canopy elements may provide direct representationof the short-circuiting processes. Overall, the spectra from both models are almost identical
for the range of the longest wavelengths down to a wavelength of approximately (2/3)h.

3.4 Two-Point Correlation Function
The mean spatial structure of coherent motions in canopy turbulence can be deduced byexamining spatial two-point correlation functions (Shaw et al. 1995; Su et al. 2000; Finnigan
and Shaw 2000), such as that between flow variables obtained at arbitrary heights and areference height, a x + r x , y + r y , z b (x, y, zr )Rab r x , r y , z, zr 
,
(31)
1/2

a 2 (z) b 2 (zr )where a and b are one of the flow variables (e.g., velocity, pressure), r x and r y are streamwise
and lateral separations, respectively, and zr is the reference height, which is fixed at a height
just above the canopy as zr /h  1.05. In the periodic domain used in the present simulations,the covariance and variances in the above equation can be calculated via the Fourier transform,

 ∗



a x + r x , y + r y , z b (x, y, zr )  F −1 b k x , k y , zr /N x N y ,(32)
a kx , k y , z 
where a and b are the Fourier transform of a and b , respectively, k x and k y are horizontal
wavenumbers, F −1 [ ] represents the inverse Fourier transform (the asterisk * denotes thecomplex conjugate), and N x N y is the total number of horizontal grid points.
As the power spectra are almost identical in the range of low wavenumbers where thespectra have significant power (Fig. 3), the two-point autocorrelation functions of the velocity
components can be foreseen as similar between the two models, as confirmed in Fig. 4. In
this figure, the lattice Boltzmann results are shown by line contours and the Navier–Stokes
results are represented by coloured shading. As expected, the two-point correlation functions
from both almost coincide. The autocorrelation of streamwise velocity component exhibitsa streamwise-elongated (Fig. 4c) and vertically tilted large-scale structure, in which the 0.1
contour reaches higher than 6 h and extends downstream beyond the figure (Fig. 4a). Suchspatial coherence in the vertical velocity component is much more compact in a horizontal
slice (Fig. 4d) and it is nearly in phase at all levels from the ground to the height 2–3 h(Fig. 4b). As these figures are similar to what has been reported following earlier studies
using wind-tunnel experiments (Shaw et al. 1995) and an LES analysis (Su et al. 2000) offlows within and above different plant canopies, the key features of the coherent structures
depicted here are considered not strongly affected by the details of canopy morphology.
Similarly to Su et al. (2000), a region of negative correlations in the function Ruu is seen
near the ground in Fig. 4a. The occurrence of such a region is probably attributable to an
adverse pressure gradient induced by large-scale sweep motions impinging on the canopy
constrained by solid ground. For more discussions, the readers should refer to Su et al. (2000)
and Shaw et al. (2013).
As suggested by Finnigan and Shaw (2000) and supported by the positive streamwise
skewness shown in Fig. 1d, sweep motions generally have greater kinetic energy than ejection motions near the canopy top. Hence, the distribution of two-point correlations, with the

123

48

T. Watanabe et al.

Fig. 3 One-dimensional power spectra at four different heights of z/h  3.05 (a, b), 2.05 (c, d), 1.05 (e, f),and 0.55 (g, h). a, c, e, and g are streamwise spectra, while b, d, f, and h are lateral spectra. Lines represent
the lattice Boltzmann results and marks denote the results of the Navier–Stokes model; L x and L y are thestreamwise and lateral dimensions of the computational domain; λx and λ y are the streamwise and lateralwavelengths

123

Large-Eddy Simulation of Neutrally-Stratified Turbulent Flow49

Fig. 3 continued

reference point near the canopy top, tends to reflect the structure of sweep motions selectively. Conducting a conditional sampling analysis on the LES results, Watanabe (2004,2009) reported that sharp streamwise gradients of the streamwise velocity component near
the canopy top, associated with a scalar micro-front between an upstream sweep and adownstream ejection, are observed most often with the passage of a streamwise-elongatedand vertically-tilted structure similar to that shown in Fig. 4. Similar spatial structures
were also detected in the LES results of Fitzmaurice et al. (2004), who used an occurrence of a positive pressure pulse near the canopy top as a trigger of conditional sampling.
As the combination of an upstream sweep and a downstream ejection causes streamwise
convergence (i.e., ∂u/∂ x < 0) near the canopy top, a positive pressure zone that develops there gives rise to divergent flow in both lateral and vertical directions. An intense
positive pressure perturbation can thus be used as an indicator of the presence of the sweep–ejection system that also causes the scalar micro-front. Therefore, although the detection
method was different between Fitzmaurice et al. (2004) and Watanabe (2004, 2009), thesweep part of the spatial patterns shown by these authors essentially represent the common
structure that induces the most energetic perturbations in the streamwise velocity component, of which the spatial structure is depicted by the two-point autocorrelation Ruu in
Fig. 4.
As shown in Fig. 4c, d, the horizontal distribution of the two-point correlation functions
exhibits lateral symmetry attributable to the statistical homogeneity in the lateral direction.
However, this picture is not quite representative because individual eddies generally haveasymmetric shapes with different orientations. To obtain averaged images of such asymmetrical spatial structure of the dominant eddies, the two-point correlation functions arerecalculated using only a part of the Fourier components in the first (k x > 0, k y > 0) andthird (k x < 0, k y < 0) quadrants in the horizontal wavenumber plane, thereby summing only
the components that have negative inclination angles of wave crest line with respect to thestreamwise direction.
Horizontal distributions of the recalculated functions Ruu , Rvu , Rwu are shown in Fig. 5,
with only the lattice Boltzmann results shown for brevity. Figure 5a reveals that the symmetricpattern seen in the streamwise velocity component (Fig. 4c) is comprised of the superposition of laterally meandering structures that are much longer streamwise than shown in
Fig. 4c. It should be noted that the unused Fourier components in the second and fourth

123

50

T. Watanabe et al.

Fig. 4 Spatial distribution of two-point autocorrelation functions of the streamwise (a, c) and vertical (b,d) velocity components in a vertical–streamwise cross-section (a, b) and a horizontal cross-section (c, d) at
the height z/h  1.05. Line contours with an interval of 0.1 denote the lattice Boltzmann results and coloured
shading represents the Navier–Stokes results123

Large-Eddy Simulation of Neutrally-Stratified Turbulent Flow

51

Fig. 4 continued

quadrants would produce mirrored patterns (with respect to the horizontal axis at r y  0) of
Fig. 5 and the summation of all components would reduce to Fig. 4. The cross-correlation
functions Rvu and Rwu (Fig. 5b, c) indicate that high values of the streamwise velocity component are accompanied not only by the negative values of the vertical velocity
component corresponding to sweep motions but also by excursions of the lateral velocity component. The lateral motions are directed systematically towards the region of low
values of the streamwise velocity component from the high-speed counterparts. Similar
results were reported by Perret and Ruiz (2013), who conducted conditional sampling of
sweep and ejection motions from wind-tunnel particle-image-velocimetry data using a biased
threshold towards positive values of the lateral velocity component. Moreover, comparedwith the streamwise-elongated feature of the streamwise velocity component, the lateraland vertical velocity components generally have shorter oblique structures that concentratenear the meandering core of the streamwise velocity component. Thus, the most energetic
sweep event, in the meandering field of streamwise velocity component, corresponds to theoblique front where a fast-moving high-speed perturbation overtakes and pushes aside lowspeed perturbations that move more slowly; this finding is consistent with Watanabe (2009).

3.5 Spectral Energy Budget of the Streamwise Structure
Although the mechanism leading to the occurrence of the streaky structures mentioned above
remains a matter of debate (Raupach et al. 1996; Finnigan et al. 2009; Watanabe 2009; Shawet al. 2013) and beyond the scope here, some insights into the process that maintains such
structures can be obtained by examining the TKE budget in more detail. Figure 6a shows the
pre-multiplied spectral distribution of the terms in the budget equation for the variance u 2 /2
calculated just above the canopy (z  1.1h). The pre-multiplied versions of the streamwise

123

52

T. Watanabe et al.

Fig. 5 Horizontal distribution at z/h  1.05 of the two-point correlation functions a Ruu , b Rvu , and c Rwu
recalculated using a limited number of Fourier components. Only the lattice Boltzmann results are shown andthe contour interval is 0.1

123Large-Eddy Simulation of Neutrally-Stratified Turbulent Flow

53

Fig. 6 Streamwise one-dimensional spectral distribution. a Pre-multiplied distribution of individual terms inthe budget equation of u u /2. Ps : shear production, Tt : turbulent transport, Rp : pressure–strain-rate correlation, TSGS : SGS transport, DSGS : SGS production. b Pre-multiplied power spectra of the u and w velocity
components. c Vertical coherence length for the u and w velocity components. Data shown in (a) and (b) areevaluated just above the canopy (z/h  1.1). The horizontal axis is the streamwise wavelength normalized by
the canopy height. Only the lattice Boltzmann results are shown

one-dimensional power spectra of the u and w components are also shown in Fig. 6b forcomparison.
As has been widely accepted, the kinetic energy is first injected into the streamwise
energy through shear production, i.e., streamwise fluctuation is induced by the vertically

123

54

T. Watanabe et al.

advected (displaced) mean velocity. The peak of shear production −Re u∗w
 du/dz (Redenotes the real part) is located in between the peaks of the pre-multiplied power spectraof the u and w components. Thus, the fundamental source of TKE includes both longand short-wavelength components. In fact, the long components (wavelengths > 10h) areresponsible for 46% of the total shear stress u w at the canopy top. However, the kinetic
energy induced by shortwave components is redistributed largely to the lateral and verticalcomponents of TKE through the action of pressure, as indicated by the pressure–strain-ratecorrelation term Im k x pu ∗ (Im denotes the imaginary part) in Fig. 6a. This is attributable
to the sharp streamwise gradients in the shortwave fluctuations of the streamwise velocity
component, which induce intense pressure perturbations to satisfy the constraint of flow
continuity. One may recall that the pressure-strain-rate correlation is the only substantial
source in the budget of lateral and vertical components of TKE, and the sum of the terms in
all three components almost vanishes because of the flow continuity. As a result, the power
spectra of the v and w components have significant values only in the short-wavelengthor large-wavenumber range (Figs. 3 and 6b), while the longwave energy tends to remainin the streamwise velocity component. The shear production to a longwave component isbalanced by the turbulent transport, which is composed of the vertical transport towards the
canopy layer and the energy cascade to shortwave components, and, to a lesser extent, direct
consumption for SGS production. Because these non-linear processes involve the higherorder interactions between the turbulent velocity fluctuations, these processes are slower
than the pressure-redistribution process, which is dominated by the linear rapid pressurethat responds instantaneously to the shear production (Pope 2000). Consequently, longwave
structures in the streamwise velocity component persist for longer.
The spectral distribution of the vertical coherence length (Fig. 6c) indicates thatstreamwise-elongated structures are generally tall. The coherence length, which is a spectral
analogue of the integral length scale, is calculated as
 z top
L αα

(33)
,
z
Cαα (k x , z, zr )dz,
(k
)
x
r
z
h

where Cαα is the coherence function between individual Fourier modes obtained at arbitraryheights and the reference height, i.e.,∼Cαα (k x , z, zr )  ∼∼∗u α (k x , y, z)u α (k x , y, zr )

2 ∼ 
2
u α (k x , y, z) u α k x , y, z r∼

1/2 .

(34)

Here, u α is the streamwise one-dimensional Fourier transform of u α , and repeated indicesin Eqs. 33 and 34 do not imply summation. The reference height is set to zr /h  1.05, andthe integral in Eq. 33 is calculated only in the region above the canopy. It should be notedthat the above definition is based on the “coherence height” used by Jiménez et al. (2004) but
differs slightly from the original definition. Figure 6c clearly indicates that the vertical lengthscales for both the streamwise and vertical velocity components increase with the streamwise
wavelength λx in the shorter range (λx ≤ 10h), where the structure of the vertical velocity
component is generally taller than that of the streamwise velocity component. For longer
wavelengths (λx > 10h), however, the length scale of the vertical velocity component levels
off, whereas the length scale of the streamwise velocity component continues to grow to thedimension of 2–3 h at the longest wavelengths. Thus, the longest modes of the streamwise123

Large-Eddy Simulation of Neutrally-Stratified Turbulent Flow55

velocity component are coherent across the roughness sublayer, which is considered to extend
to a height of roughly 3h.
Taken altogether, the larger contribution of longwave components to the streamwise
velocity component (Fig. 6b) results in a more streamwise-elongated, taller coherent structure (Fig. 4a) in comparison with the compact structure of the vertical velocity component(Fig. 4b).

4 Conclusions
We have described a central-moments-based lattice Boltzmann method for the simulationof neutrally-stratified turbulent flows with the LES approach. By conducting comparative
simulations of the airflow within and above a homogeneous plant canopy, the performanceof the model was evaluated against a conventional LES model based on the Navier–Stokes
equations. The simulated turbulence statistics such as mean velocity, velocity variances,velocity skewness, and power spectra are almost identical between the two models. The
spatial structure of the coherent eddies and their maintenance processes were confirmed
properly represented by the lattice Boltzmann approach through analysis of the TKE budget
and the spatial distribution of two-point correlation functions.
As already shown by previous studies, the most prominent spatial structures in thestreamwise velocity component present a streaky form in the roughness sublayer over bothvegetation and urban canopies. An examination of the spectral TKE budget shows that thedifference in the pressure-redistribution process between short- and long-wavelength components of the shear-generated streamwise energy plays a role in the maintenance of thestreaky structures. Although the turbulent pressure rapidly redirects short-wavelength components of streamwise energy to the lateral and vertical velocity components, long-wavelengthcomponents of energy tend to remain in the streamwise velocity component because thepressure-redistribution process works only marginally because of the small streamwise gradients in the long-wavelength components.
Based on our experience, the lattice Boltzmann algorithm is quite simple and its coding isstraightforward. Its particulate nature and local dynamics are suitable for high-performancecomputations of the atmospheric-boundary-layer flows. Although we only dealt with a neutrally stratified flow, the method can be applied to non-neutral conditions by including thebuoyancy in the forcing terms, with the aid of a coupled solver for (potential) temperature. The
temperature solver could also use a lattice Boltzmann approach for predicting the distributionfunction of the internal energy (e.g., He et al. 1998a; Wang et al. 2018) or a conventional
finite-difference approach for the macroscopic thermodynamic equation.
Acknowledgements This study was supported in part by the Grant for Joint Research Program of the Institute of Low Temperature Science, Hokkaido University (18G009, 19G007). The reference simulation was
performed on the supercomputing system (NEC SX-9) at the Computer Centre for Agriculture, Forestry and
Fisheries Research Information Technology Centre, MAFF Japan during a joint research project supported by
JSPS KAKENHI Grant Number JP19510020. Thanks are also due to Dr. Naoyuki Onodera for his invaluable
discussions.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give
appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence,and indicate if changes were made. The images or other third party material in this article are included in thearticle’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is
not included in the article’s Creative Commons licence and your intended use is not permitted by statutory

123

56

T. Watanabe et al.

regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix: Equations for the Collision and Forcing
The post-collision central moments are generally calculated by Eq. 2 with the collision
process parametrized as a relaxation towards equilibrium (Eq. 13). However, to maintainthe rotational invariance, some moments need to be rearranged before relaxation, and somerelaxation coefficients must be assigned the same value (Geier et al. 2015). The actually
implemented equations are described here.
For the conservative moments (i.e., the zeroth- and first-order moments, related to themass and momentum conservations, respectively), the relaxation coefficients must be zero,and the post-collision moments are given as∗κ000 κ000 ,

(35)

∗
κ100
 κ100 + Fx  −κ100 ,

(36a)

∗
 κ010 + Fy  −κ010 ,
κ010

(36b)

∗
 κ001 + Fz  −κ 001 .
κ001

(36c)

The second-order moments are rearranged to correctly recover the traceless and trace components of the momentum flux tensor, to which the coefficients ω1 and ω2 are applied
respectively,∗
κ110
 (1 − ω1 )κ110 ,

(37a)

∗
 (1 − ω1 )κ011 ,
κ011

(37b)

∗
 (1 − ω1 )κ101 ,
κ101



ω1  2 ∂u
∂v
∗
∗
κ200
u
− v2
,
− κ020 (1 − ω1 )(κ200 − κ020 ) − 3ρ 1 −2  ∂x
∂y 


ω1
∂u
∂w
∗
∗
− κ002
 (1 − ω1 )(κ200 − κ002 ) − 3ρ 1 −
κ200
u2
− w2
,
2
∂x
∂z
∗
∗
∗
κ200
+ κ020+ κ002
 (1 − ω2 ) (κ200 + κ020 + κ002 ) + ω2 ρ



∂v
∂w
ω2  2 ∂uu.
− 3ρ 1 −
+ v2
+ w2
2
∂x
∂y
∂z

(37c)
(38a)
(38b)

(39)

Final terms in Eqs. 38 and 39 are partial corrections for the lack of Galilean invariance ofthe original lattice Boltzmann equation (Geier et al. 2015), and the velocity gradients usedin these terms are calculated by Eq. 48 shown below. Similarly, the post-collision values of
the higher-order moments are calculated by the rearranged equations asω3  2Fx∗∗
κ120+ κ102
 (1 − ω3 )(κ120 + κ102 ) + 1 −
,
(40a)
2
3


ω3 2Fy∗∗,(40b)+ κ012
 (1 − ω3 )(κ210 + κ012 ) + 1 −
κ210
2
3


ω3 2Fz∗∗+ κ021 (1 − ω3 )(κ201 + κ021 ) + 1 −
κ201
,
(40c)
2
3

123

Large-Eddy Simulation of Neutrally-Stratified Turbulent Flow

57

∗
∗
κ120
− κ102
 (1 − ω4 )(κ120 − κ102 ),

(41a)

∗
∗
κ210
− κ012
 (1 − ω4 )(κ210 − κ012 ),

(41b)

∗
∗
κ201
− κ021 (1 − ω4 )(κ201 − κ021 ),

(41c)

∗
κ111
 (1 − ω5 )κ111 ,

(42)

∗
∗
∗
κ220
− 2κ202+ κ022
 (1 − ω6 )(κ220 − 2κ202 + κ022 ),

(43a)

∗
∗
∗
κ220
+ κ202− 2κ022
 (1 − ω6 )(κ220 + κ202 − 2κ022 ),

(43b)

ρ
∗
∗
∗
κ220
+ κ202+ κ022
 (1 − ω7 )(κ220 + κ202 + κ022 ) + ω7 ,
3

(44)

∗
κ211
 (1 − ω8 )κ211 ,

(45a)

∗
κ121
 (1 − ω8 )κ121 ,

(45b)

∗
κ112
 (1 − ω8 )κ112 ,
ω9  Fz∗ (1 − ω9 )κ221 + 1 −
κ221
,
2 9


ω9 Fy
∗
 (1 − ω9 )κ212 + 1 −
κ212
,
2 9


ω9 Fx
∗
,
 (1 − ω9 )κ122 + 1 −
κ122
2 9
ρ
∗
κ222
 (1 − ω10 )κ222 + ω10 ,
27

(45c)
(46a)
(46b)
(46c)
(47)

where ω3,4,5 are the relaxation coefficients for the third-order moments, ω6,7,8 are those forthe fourth-order moments, and ω9 and ω10 are those for the fifth- and the sixth-order moments,
respectively. The forcing terms in the higher-order moments are included in Eqs. 40 and 46,
as suggested by Fei and Luo (2017) and Fei et al. (2018). The velocity gradients in Eqs. 38
and 39 are evaluated from the second-order central moments as (Geier et al. 2017)∂uω2ω1
 − (2κ200 − κ020 − κ002 ) −(κ200 + κ020 + κ002 − ρ),∂x2ρ2ρ∂u 3ω1∂v
+(κ200 − κ020 ),
∂y
∂x
2ρ
∂w
∂u 3ω1+(κ200 − κ002 ).
∂z
∂x
2ρ

(48a)
(48b)
(48c)

Other components of the strain rates are similarly evaluated as
∂v ∂u
3ω1
+
−
κ110 ,
∂x ∂y
ρ
∂w ∂u
3ω1
+
−
κ101 ,
∂x
∂z
ρ
3ω1
∂w ∂v
+
−
κ011 .
∂ y ∂z
ρ

(49a)
(49b)
(49c)

The final three relationships are used in the evaluation of the subgrid eddy viscosity (Eq. 15).

123

58

T. Watanabe et al.

References
Ahmad NH, Inagaki A, Kanda M, Onodera N, Aoki T (2017) Large-eddy simulation of the gust index in an
urban area using the lattice Boltzmann method. Boundary-Layer Meteorol 163:447–467
Aidum CK, Clausen JR (2010) Lattice-Boltzmann method for complex flows. Annu Rev Fluid Mech
42:439–472
Bergstöm H, Högström U (1989) Turbulent exchange above a pine forest. II Organized structures. BoundaryLayer Meteorol 49:231–263
Bohrer G, Katul GG, Walko RL, Avissar R (2009) Exploring the effects of microscale structural heterogeneityof forest canopies using large-eddy simulations. Boundary-Layer Meteorol 132:351–382
Cava D, Katul GG (2008) Spectral short-circuiting and wake production within the canopy trunk space of an
Alpine hardwood forest. Boundary-Layer Meteorol 126:415–431
Chen S, Doolen GD (1998) Lattice Boltzmann method for fluid flows. Annu Rev Fluid Mech 30:329–364
Coceal O, Thomas TG, Catro IP, Belcher SE (2006) Mean flow and turbulence statistics over groups ofurban-like cubical obstacles. Boundary-Layer Meteorol 121:491–519
Delbosc N, Summers JL, Khan AI, Kapur N, Noakes CJ (2014) Optimized implementation of the lattice
Boltzmann method on a graphics processing unit towards real-time fluid simulation. Comput Math Appl
67:462–475
Denmead OT, Bradley EF (1985) Flux-gradient relationships in a forest canopy. In: Hutchison BA, Hicks BB(eds) The forest–atmosphere interaction. Springer, Dordrecht, pp 421–442
Dupont S, Brunet Y (2009) Coherent structures in canopy edge flow: a large-eddy simulation study. J Fluid
Mech 630:93–128
Dupont S, Irvine MR, Bonnefond JM, Lamaud E, Brunet Y (2012) Turbulent structures in a pine forest witha deep and sparse trunk space: stand and edge regions. Boundary-Layer Meteorol 143:309–336
Dwyer MJ, Patton EG, Shaw RH (1997) Turbulent kinetic energy budgets from a large-eddy simulation of
airflow above and within a forest canopy. Boundary-Layer Meteorol 84:23–43
Fei L, Luo KH (2017) Consistent forcing scheme in the cascaded lattice Boltzmann method. Phys Rev E
96:053307
Fei L, Luo KH, Li Q (2018) Three-dimensional cascaded lattice Boltzmann method: improved implementationand consistent forcing scheme. Phys Rev E 97:053309
Finnigan J (2000) Turbulence in plant canopies. Annu Rev Fluid Mech 32:519–571
Finnigan JJ, Shaw RH (2000) A wind-tunnel study of air flow in waving wheat: an EOF analysis of the structureof the large-eddy motion. Boundary-Layer Meteorol 96:211–255
Finnigan JJ, Shaw RH, Patton EG (2009) Turbulence structure above a vegetation canopy. J Fluid Mech
637:387–424
Fitzmaurice L, Shaw RH, Paw UKT, Patton EG (2004) Three-dimensional scalar microfront systems in alarge-eddy simulation of vegetation canopy flow. Boundary-Layer Meteorol 112:107–127
Gao W, Shaw RH, Paw UKT (1989) Observation of organized structure in turbulent flow within and above aforest canopy. Boundary-Layer Meteorol 47:349–377
Geier M, Greiner A, Korvink JG (2006) Cascaded digital lattice Boltzmann automata for high Reynolds
number flow. Phys Rev E 73:066705
Geier M, Pasquali A, Schönherr M (2017) Parametrization of the cumulant lattice Boltzmann method for
fourth order accurate diffusion. Part I: derivation and validation. J Comput Phys 348:862–888
Geier M, Schönherr M, Pasquali A, Krafczyk M (2015) The cumulant lattice Boltzmann equation in three
dimensions: theory and validation. Comput Math Appl 70:507–547
Guo Z, Zhao TS (2002) Lattice Boltzmann model for incompressible flows through porous media. Phys Rev
E 66:036304
Guo Z, Zheng C, Shi B (2002) Discrete lattice effects on the forcing term in the lattice Boltzmann method.
Phys Rev E 65:046308
He X, Chen S, Doolen GD (1998a) A novel thermal model for the lattice Boltzmann method in incompressible
limit. J Comput Phys 146:282–300
He X, Shan X, Doolen GD (1998b) Discrete Boltzmann equation model for nonideal gases. Phys Rev E
57:R13–R16
Inagaki A, Casillo MCL, Yamashita Y, Kanda M, Takimoto H (2012) Large-eddy simulation of coherent flow
structures within a cubical canopy. Boundary-Layer Meteorol 142:207–222
Inagaki A, Kanda M, Ahmad NH, Yagi A, Onodera N, Aoki T (2017) A numerical study of turbulence statistics
and the structure of a spatially-developing boundary layer over a realistic urban geometry. BoundaryLayer Meteorol 164:161–181
Inagaki M, Nagaoka M, Horinouchi N, Suga K (2010) Large eddy simulation analysis of engine steady intake
flows using a mixed-time-scale subgrid-scale model. Int J Engine Res 11:229–241

123

Large-Eddy Simulation of Neutrally-Stratified Turbulent Flow

59

Jiménez J, del Álamo JC, Flores O (2004) The large-scale dynamics of near-wall turbulence. J Fluid Mech
505:179–199
Kaimal JC, Finnigan JJ (1994) Atmospheric boundary layer flows. Oxford University Press, New York
Kanda M (2006) Large-eddy simulations on the effects of surface geometry of building arrays on turbulent
organized structures. Boundary-Layer Meteorol 118:151–168
Kanda M, Hino M (1994) Organized structures in developing turbulent flow within and above a plant canopy,
using a large eddy simulation. Boundary-Layer Meteorol 68:237–257
Kanda M, Moriwaki R, Kasamatsu F (2004) Large-eddy simulation of turbulent organized structures within
and above explicitly resolved cube arrays. Boundary-Layer Meteorol 112:343–368
Katul GG, Cava D, Siqueira M, Poggi D (2013) Scalar turbulence within the canopy sublayer. In: Venditti JG,
Best JL, Church M, Hardy RJ (eds) Coherent flow structures at Earth’s surface. Wiley, Oxford, pp 73–95
Kobayashi H (2005) The subgrid-scale models based on coherent structures for rotating homogeneous turbulence and turbulent channel flow. Phys Fluids 17:045104
Krüger T, Kusumaatmaja H, Kuzmin A, Shardt O, Silva G, Viggen EM (2017) The lattice Boltzmann method.
Springer, Bern
Lallemand P, Luo LS (2003) Lattice Boltzmann method for moving boundaries. J Comput Phys 184:406–421
Meyers TP, Baldocchi DD (1991) The budgets of turbulent kinetic energy and Reynolds stress within andabove a deciduous forest. Agric For Meteorol 53:207–222
Mohamad AA (2011) Lattice Boltzmann method. Springer, London
Onodera N, Aoki T, Shimokawabe T, Kobayashi H (2013) Large-scale LES wind simulation using lattice
Boltzmann method for a 10 km × 10 km area in metropolitan Tokyo. Tsubame ESJ 9:2–8
Ouwersloot HG, Moene AF, Attema JJ, De Arellano JVG (2017) Large-eddy simulation comparison of neutral
flow over a canopy: sensitivities to physical and numerical conditions, and similarity to other representations. Boundary-Layer Meteorol 162:71–89
Patton EG, Shaw RH, Judd MJ, Raupach MR (1998) Large-eddy simulation of windbreak flow. BoundaryLayer Meteorol 87:275–307
Patton EG, Sullivan PP, Shaw RH, Finnigan JJ, Weil JC (2016) Atmospheric stability influences on coupledboundary layer and canopy turbulence. J Atmos Sci 73:1621–1647
Perret L, Ruiz T (2013) SPIV analysis of coherent structures in a vegetation canopy model flow. In: Venditti
JG, Best JL, Church M, Hardy RJ (eds) Coherent flow structures at Earth’s surface. Wiley, Oxford, pp
161–174
Poggi D, Porporato A, Ridolfi L, Albertson JD, Katul GG (2004) The effect of vegetation density on canopy
sub-layer turbulence. Boundary-Layer Meteorol 111:565–587
Pope SB (2000) Turbulent flows. Cambridge Univ Press, Cambridge
Premnath KN, Banerjee S (2009) Incorporating forcing terms in cascaded lattice Boltzmann approach bymethod of central moments. Phys Rev E 80:036702
Premnath KN, Banerjee S (2011) On the three-dimensional central moment lattice Boltzmann method. J Stat
Phys 143:747–794
Raupach MR, Coppin PA, Legg BJ (1986) Experiments on scalar dispersion within a model plant canopy part
I: the turbulence structure. Boundary-Layer Meteorol 35:21–52
Raupach MR, Finnigan JJ, Brunet Y (1996) Coherent eddies and turbulence in vegetation canopies: themixing-layer analogy. Boundary-Layer Meteorol 78:351–382
Raupach MR, Shaw RH (1982) Averaging procedures for flow within vegetation canopies. Boundary-Layer
Meteorol 22:79–90
Raupach MR, Thom AS (1981) Turbulence in and above plant canopies. Annu Rev Fluid Mech 13:97–129
Rinaldi PR, Dari EA, Vénere MJ, Clausse A (2012) A lattice-Boltzmann solver for 3D fluid simulation on
GPU. Simul Model Pract Theory 25:163–171
Seginer I, Mulhearn PJ, Bradley EF, Finnigan JJ (1976) Turbulent flow in a model plant canopy. BoundaryLayer Meteorol 10:423–453
Shaw RH (1977) Secondary wind speed maxima inside plant canopy. J Appl Meteorol 16:514–521
Shaw RH, Brunet Y, Finnigan JJ, Raupach MR (1995) A wind tunnel study of air flow in waving wheat:two-point velocity statistics. Boundary-Layer Meteorol 76:349–376
Shaw RH, Patton EG, Finnigan JJ (2013) Coherent eddy structures over plant canopies. In: Venditti JG, BestJL, Church M, Hardy RJ (eds) Coherent flow structures at Earth’s surface. Wiley, Oxford, pp 149–159
Shaw RH, Schumann U (1992) Large-eddy simulation of turbulent flow above and within a forest. BoundaryLayer Meteorol 61:47–64
Shen S, Leclerc MY (1997) Modelling the turbulence structure in the canopy layer. Agric Forest Meteorol
87:3–25
Su HB, Shaw RH, Paw UKT (2000) Two-point correlation analysis of neutrally stratified flow within andabove a forest from large-eddy simulation. Boundary-Layer Meteorol 94:423–460

123

60

T. Watanabe et al.

Su HB, Shaw RH, Paw UKT, Moeng CH, Sullivan PP (1998) Turbulent statistics of neutrally stratified flow
within and above a sparse forest from large-eddy simulation and field observations. Boundary-Layer
Meteorol 88:363–397
Suga K, Kuwata Y, Takashima K, Chikasue R (2015) A D3Q27 multiple-relaxation-time lattice Boltzmann
method for turbulent flows. Comput Math Appl 69:518–529
Wang Y, MacCall BT, Hocut CM, Zeng X, Fernando HJS (2018) Simulation of stratified flows over a ridge
using a lattice Boltzmann model. Environ Fluid Mech. https://doi.org/10.1007/s10652-018-9599-3
Watanabe T (1993) The bulk transfer coefficients over a vegetated surface based on K-theory and a 2nd-orderclosure model. J Meterol Soc Jpn 71:33–42
Watanabe T (2004) Large-eddy simulation of coherent turbulence structures associated with scalar ramps overplant canopies. Boundary-Layer Meteorol 112:307–341
Watanabe T (2009) LES study on the structure of coherent eddies inducing predominant perturbations invelocities in the roughness sublayer over plant canopies. J Meteorol Soc Jpn 87:39–56
Wilson NR, Shaw RH (1977) A higher order closure model for canopy flow. J Appl Meteorol 16:1197–1205
Yue W, Parlange MB, Moneveau C, Zhu W, Van Haout R, Katz J (2007) Large-eddy simulation of plant canopy
flows using plant-scale representation. Boundary-Layer Meteorol 124:183–203
Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps andinstitutional affiliations.

123