Boundary-Layer Meteorology (2021) 181:39–71
https://doi.org/10.1007/s10546-021-00633-1
RESEARCH ARTICLE

Coherent Eddies Transporting Passive Scalars Throughthe Plant Canopy Revealed by Large‑Eddy Simulations Using
the Lattice Boltzmann Method
Tsutomu Watanabe1 · Marie Takagi2 · Kou Shimoyama1 · Masayuki Kawashima1 ·Naoyuki Onodera3 · Atsushi Inagaki4
Received: 19 November 2020 / Accepted: 27 May 2021 / Published online: 9 July 2021
© The Author(s) 2021

Abstract
A double-distribution-function lattice Boltzmann model for large-eddy simulations of apassive scalar field in a neutrally stratified turbulent flow is described. In simulations of thescalar turbulence within and above a homogeneous plant canopy, the model’s performanceis found to be comparable with that of a conventional large-eddy simulation model basedon the Navier–Stokes equations and a scalar advection–diffusion equation in terms of themean turbulence statistics, budgets of the second moments, power spectra, and spatial twopoint correlation functions. For a top-down scalar, for which the plant canopy serves as a
distributed sink, the variance and flux of the scalar near the canopy top are predominantly
determined by sweep motions originating far above the canopy. These sweep motions,which have spatial scales much larger than the canopy height, penetrate deep inside thecanopy and cause scalar sweep events near the canopy floor. By contrast, scalar ejection
events near the canopy floor are induced by coherent eddies generated near the canopy top.
The generation of such eddies is triggered by the downward approach of massive sweep
motions to existing wide regions of weak ejective motions from inside to above the canopy.
The non-local transport of scalars from above the canopy to the canopy floor, and viceversa, is driven by these eddies of different origins. Such non-local transport has significant
implications for the scalar variance and flux budgets within and above the canopy, as wellas the transport of scalars emitted from the underlying soils to the atmosphere.
Keywords Double distribution function · Gram–Charlier distribution · Non-local
transport · Scalar turbulence · Sweep and ejection* 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

Graduate School of Environmental Science, Hokkaido University, Kita 10 Nishi 5, Kita‑ku,
Sapporo 060‑0810, Japan

3

Centre for Computational Science and E‑Systems, Japan Atomic Energy Agency, 178‑4‑4
Wakashiba, Kashiwa 277‑0871, Japan

4

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

13

Vol.:(0123456789)

40

T. Watanabe et al.

1 Introduction
A plant canopy produces a peculiar internal environment through physical and biological processes. The resulting microclimate exerts an essential influence on the energybalance, plant physiology, hydrological processes (including snowpack and frozen soil),and ecological processes in the local plant community, thereby affecting the spatio–temporal variations in exchanges of momentum and scalars such as heat, water vapour, carbon dioxide, other gasses, and particulate matter between the canopy elements and air
within the canopy layer. Turbulent eddies transporting momentum and scalars to/fromthe atmosphere above serve to renew the microclimate inside the canopy and convey thephysical and biological influences of vegetated surfaces to the atmosphere, thus bridging from the microscale interactions between plants and the microclimate to the mesoscale or larger-scale land–atmosphere interactions.
Over recent decades, numerous studies based on experiments, theoretical considerations, and numerical analyses have investigated the turbulent transport process within
and above plant canopies, and the common view that coherent eddies play a significant role in the transport of momentum and scalars has been established (Raupach andThom 1981; Finnigan 2000; Brunet 2020). The coherent eddies, which have verticalscales comparable to or often larger than the canopy height, bring about the direct transport between the canopy layer and the atmosphere above through sweep and ejection
motions. Such non-local transport phenomena, the occurrence of which is identified
by vertically coherent ramp signals in scalar time traces measured within and above
plant canopies (Bergstöm and Högström 1989; Gao et al. 1989), occur irrespective ofthe local vertical gradients of mean quantities (i.e. wind speed, scalar concentrations),
thereby causing secondary wind-speed maxima (Shaw 1977) and counter-gradient scalar fluxes (Denmead and Bradley 1985). Accordingly, the gradient–diffusion model (orK-theory) inside the canopy has been found to be inadequate (Katul et al. 2013). The
non-local transport should be treated explicitly to ensure a realistic representation of
the canopy microclimate, although the gradient–diffusion model (including the lowerlevel second-order closure models by Mellor and Yamada 1974) is adequate in scenarios
where an approximate representation of the net vertical transport suffices. These include
the bulk parametrization of momentum and scalar exchanges over vegetated surfaces(Kondo and Watanabe 1992; Watanabe 1994; Takata et al. 2003), analysis of the energybalance during rainfall interception events (Watanabe and Mizutani 1996), and modelling the dynamic interactions between the microclimate and plant community (Watanabe et al. 2004; Toda et al. 2009).
Higher-order closure models (Wilson and Shaw 1977; Meyers and Paw U 1986, 1987;
Watanabe 1993; Katul and Albertson 1998; Tanaka 2001) and Lagrangian approaches(Raupach 1989; Siqueira et al. 2000; Ueyama et al. 2014) have been adopted as alternatives to the gradient–diffusion model. In principle, the terms of third-order (or triple)
moments and the pressure transport in the higher-order closure models, as well as the‘far-field’ diffusion term in the Lagrangian approaches, represent the non-local transport. However, these terms are parametrized heuristically based on dimensional considerations and mathematical consistency and are reported to be insufficient in some situations (Katul and Chang 1999). As a thorough and precise understanding of the natureand mechanism of the non-local transport phenomena is still lacking, it is important toobtain deeper insights through analyses of the phenomenology of canopy turbulence,
as in previous studies (Raupach 1981; Katul et al. 1997a, b, 2018; Cava et al. 2006) in

13

Coherent Eddies Transporting Passive Scalars Through the Plant…

41

which the third-order moments are linked to the sweep–ejection cycle that dominates
the momentum and scalar transport (Finnigan 2000).
Although many studies have been conducted to elucidate the turbulence structure
observed near the canopy top (Raupach et al. 1996; Finnigan 2000; Brunet 2020), lessattention has been paid to the eddies that penetrate through the canopy and transport
momentum and scalars directly to/from near the canopy floor from/to above the canopy,
thereby promoting the deepest non-local transport within the canopy. A better understanding of such canopy-penetrating transport would also have important agricultural and biogeochemical implications, such as accurate predictions of water/soil temperatures underneathcrop canopies (Kuwagata et al. 2008), the specification of soil respiratory flux above a forest ecosystem using a measurement of stable isotopes (Murayama et al. 2010), emissionsof methane from permafrost covered by forests (Iwata et al. 2015), and pathways whereby
volatile organic matter emitted from forest soils influence the activities of cloud condensation nuclei in the atmosphere (Müller et al. 2017). Therefore, the present study focuses on
the nature of the coherent eddies that transport passive scalars to/from near the canopy
floor from/to the atmosphere above the canopy, as well as the eddies that cause the predominant fluctuations in the scalar concentrations near the canopy top.
As a tool for numerical simulations of the scalar turbulence within and above the canopy, a new code for the double-distribution-function lattice Boltzmann method is introduced. The lattice Boltzmann method is a relatively recent technique in computational
fluid dynamics (Chen and Doolen 1998; Aidum and Clausen 2010). The method predicts
the spatio–temporal variation of the distribution function of the microscopic fluid-particle
velocities based on the discretized Boltzmann equation, which includes the key processesof collision, forcing, and streaming. The macroscopic variables (fluid density, velocity, andpressure) derived from low-order moments of the distribution function are mathematicallyshown to satisfy the weakly compressible Navier–Stokes equations (Krüger et al. 2017).
Recently, the method has been applied to the atmospheric boundary-layer flow over urban
canopies by Onodera et al. (2013, 2021), Ahmad et al. (2017), and Inagaki et al. (2017).
For the airflow within and above a plant canopy, Watanabe et al. (2020) reported that alarge-eddy simulation (LES) using a central-moment-based lattice Boltzmann method
reproduces the turbulent velocity field as well as a conventional LES model based on the
Navier–Stokes equations. This motivated us to extend the lattice Boltzmann method to theadvection and diffusion processes of a passive scalar within and above the plant canopy.
The scalar field in a turbulent flow is expressed by a double-distribution-function model, in
which the distribution function of the scalar density is treated separately to the velocity distribution function (Shan 1997; He et al. 1998). Regarding memory consumption, the lattice
Boltzmann method uses a number of prognostic variables (i.e. 7–27 discretized distribution functions depending on the lattice system adopted) to predict a single scalar quantity.
In this respect, a conventional model based on the advection–diffusion equation for themacroscopic scalar density is more advantageous. However, the lattice Boltzmann method
employs an almost local algorithm that is suitable for parallel computations, uses a relatively simple implementation of boundary conditions, and provides an exact representationof the streaming, which enables advection calculations with higher numerical accuracy.
Therefore, the performance of this method is worth testing.
This paper describes a lattice Boltzmann model for simulating a passive-scalar field ina neutrally stratified turbulent flow within and above a plant canopy (Sect. 2). The model
is first evaluated with respect to a conventional LES model based on the Navier–Stokes
equations and the advection–diffusion equation for a scalar field (Sect. 3.1). Then, the characteristics of the eddy motions that contribute to the scalar variance and fluxes near the

13

42

T. Watanabe et al.

canopy top (Sect. 3.2) and the spatio–temporal structure of coherent eddies that promote
non-local scalar transport to/from deep inside the canopy from/to the atmosphere above are
examined (Sect. 3.3).

2 Models
A Cartesian coordinate system is adopted in which x , y, and z denote the streamwise, lateral, and vertical directions, respectively. The velocity components in these coordinates are
represented as u, v, and w. Vector and tensor notation, such as 𝐱 = ( x1, x2, x3) = (x , y, z ) and𝐮 = (u1, u2, u3) = (u, v, w), is also used. Unless otherwise stated, repeated indices in the tensor notation imply a summation over all possible values of the index.

2.1 Navier–Stokes Model
A reference simulation is performed using a LES model (Watanabe 2004, 2009), which
is based on the Navier–Stokes equations, the incompressible continuity equation, andthe advection–diffusion equation for a scalar. The subgrid-scale (SGS) turbulence kinetic
energy (TKE) is modelled using a separate prognostic equation. The subgrid eddy viscosity
is calculated from the SGS TKE, and a constant subgrid turbulent Schmidt number (= 1∕3)
(Deardorff 1971) is used to parametrize the subgrid eddy diffusivity of the scalar. The plant
canopy is represented by distributed sinks/sources of momentum and the scalar. The set-up
for the reference simulation is the same as that of the lattice Boltzmann model (describedlater), except that the exact value of the flow-driving force and the domain size are different(Watanabe et al. 2020). The velocity statistics obtained from the reference simulation werereported in our previous study (Watanabe et al. 2020), and the scalar statistics are included
herein. This reference model, including the advection–diffusion equation for a scalar, ishereafter referred to as the Navier–Stokes model.

2.2 Lattice Boltzmann Method
In the following, the equations of the lattice Boltzmann method are presented using the lattice-unit system, in which all variables are normalized by the lattice spacing Δx , timestep
Δt , reference fluid density 𝜌0, and reference scalar density s0. Therefore, in the normalized
equations, these reference scales read Δx = Δt = 𝜌0 = s0 = 1 and are eliminated whenever
there is no possibility of confusion.

2.2.1 Central‑Moment‑based Model for Airflow
The model used to simulate airflow has been fully described by Watanabe et al. (2020); asummary is repeated here. The basic framework of the model is the central-moment-basedlattice Boltzmann method in a three-dimensional 27-velocity (D3Q27) lattice system. Hereand hereafter, the notation DdQq represents the lattice type, with d indicating the number of spatial dimensions and q the number of discrete microscopic velocities considered.
The model describes the fate of the distribution function fijk (𝐱, t) of discrete microscopic
velocities (ic, jc, kc), where the indices i, j, k ∈ {−1, 0, 1} represent the direction of the discretized velocities and c = Δx∕Δt (= 1) denotes the lattice velocity unit.

13

43

Coherent Eddies Transporting Passive Scalars Through the Plant…

The distribution function and the macroscopic flow variables are updated by the following
procedure. At every timestep, the central moments of the distribution function 𝜅pqr, are calculated by𝜅pqr (𝐱, t) =

1
1
1
∑
∑
∑

(ic − u)p (jc − v)q (kc − w)r fijk (𝐱, t),(1)

i=−1 j=−1 k=−1where the indices p, q, r ∈ {0, 1, 2} are the directional orders of moments, and the sum of
the indices ( p + q + r ) defines the formal order of the moments. Using the multi-relaxation-time collision scheme and a forcing scheme applied in the central-moment space,∗ is obtained in terms of the pre-collisional central
the post-collisional central moment 𝜅pqr()
moments, the macroscopic fluid density and velocity, a net external force 𝐅 = Fx , Fy , F𝐳working on the fluid, and the collisional relaxation coefficients 𝜔1 , … , 𝜔10 as
)
(
∗
𝜅pqr
(𝐱, t) = func 𝛋, 𝜌, 𝐮, 𝐅, 𝜔1 , … , 𝜔10 .
(2)
The functional form of Eq. 2 is described in Watanabe et al. (2020). The inverse transform∗, which is transof Eq. 1 is then used to calculate the post-collision distribution function fijk
lated to the adjacent grid points in the next streaming step as∗fijk (x + icΔt, y + jcΔt, z + kcΔt, t + Δt) = fijk(x, y, z, t).

(3)

After the application of appropriate boundary conditions to this new distribution function,
the macroscopic variables at the next timestep are calculated from

𝜌=

1
1
1
∑
∑
∑

(4)

fijk ,

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

𝜌u =

1
1
1
∑
∑
∑

(ic)f ijk +i=−1 j=−1 k=−1

𝜌v =

1
1
1
∑
∑
∑

(jc)fijk +i=−1 j=−1 k=−1

𝜌w =

1
1
1
∑
∑
∑

(kc)fijk +i=−1 j=−1 k=−1

p=

𝜌
,
3

Fx
,
2
Fy

(5a)

,

(5b)

Fz
,
2

(5c)

2

(6)

where p is the macroscopic static pressure.
To enable simulations of high-Reynolds-number flows, the stress exerted by the SGS
motion is modelled as)(𝜕u𝛼 𝜕u𝛽2++ kSGS 𝛿𝛼𝛽 ,𝜏𝛼𝛽 = −𝜌𝜐SGS
(7)
𝜕x𝛽
𝜕x𝛼
3

13

T. Watanabe et al.

44

where 𝜏𝛼𝛽 represents the SGS stress, the indices 𝛼 and 𝛽 indicate the spatial coordinates,
𝜐SGS denotes the subgrid eddy viscosity, and 𝛿𝛼𝛽 is the Kronecker delta. The SGS kinetic
energy denoted by kSGS is estimated from()2
kSGS = 𝜌ck u𝛼 − u𝛼 ,(8)where ck (= 1) is a model constant and u𝛼 is a filtered velocity component obtained by

u𝛼 =

)
1( B
E
N
T
u + uS𝛼 + uW
𝛼 + 6u𝛼 + u𝛼 + u𝛼 + u𝛼 ,
12 𝛼

(9)

in which uB𝛼 , uS𝛼 , … denote the velocities at the six surrounding grid points (Suga et al.
2015). The subgrid eddy viscosity is written in terms of the resolved strain-rate tensor S𝛼𝛽as√𝜐SGS = CSGS (Δx)2 2S𝛼𝛽 S𝛼𝛽 ,(10)where CSGS is a parameter determined according to the coherent-structure Smagorinsky
model of Kobayashi (2005). The evaluated value of 𝜐SGS is reflected in the lattice Boltzmann algorithm through the relaxation coefficient 𝜔1 for a group of second-order central
moments (Watanabe et al. 2020) as
(
)
1 1
1
𝜐 + 𝜐SGS =−.
(11)
3 𝜔1 2The other relaxation coefficients of the central-moment-based collision model are set as𝜔2 = 𝜔6 = 𝜔7 = 𝜔8 = 𝜔10 = 𝜔1 and 𝜔3 = 𝜔4 = 𝜔5 = 𝜔9 = 1 to reduce the spurious
numerical viscosity as simply as possible (Geier et al. 2015).
With these parametrizations, the lattice Boltzmann airflow model predicts the macroscopic
density and velocity that satisfy the filtered continuity equation𝜕𝜌 𝜕𝜌u𝛼
+
=0
𝜕t
𝜕x𝛼

(12)

and the filtered Navier–Stokes equations([)]) 𝜕u𝛼 𝜕u𝛽(𝜕𝜌u𝛼 𝜕𝜌u𝛼 u𝛽
𝜕p
𝜕
2
+
=−
+
𝜌 𝜐 + 𝜐SGS+
− kSGS 𝛿𝛼𝛽 + F𝛼 . (13)𝜕t𝜕x𝛽
𝜕x𝛼 𝜕x𝛽𝜕x𝛽
𝜕x𝛼
32.2.2 Single‑Relaxation‑Time Model for a Scalar Field
The model for a scalar field is based on a single-relaxation-time model, i.e. theBhatnagar–Gross–Krook model. This is the most common and simplest collision model
used in the lattice Boltzmann literature (Chen and Doolen 1998; Aidum and Clausen2010). A D3Q19 lattice is adopted as a trade-off between numerical accuracy and computational cost—considerable shortwave noise appeared in the power spectra of the scalardensity in preliminary simulations using more efficient D3Q7 and D3Q15 lattices, while
the most accurate D3Q27 lattice requires more memory resources.

13

Coherent Eddies Transporting Passive Scalars Through the Plant…

45

In this framework, the spatio–temporal development of the distribution function of ascalar, gl (l = 0, … , 18), is described by (Krüger et al. 2017)
(
)
(
)
1
1
eq
gl (𝐱, t) + Ql (𝐱, t) − gl (𝐱, t) + Ql (𝐱, t), (14)gl 𝐱 + 𝐜l Δt, t + Δt = gl (𝐱, t) −𝜏2gwhere 𝐜l is one of the 19 microscopic velocities, defined as 𝐜0 = (0, 0, 0), 𝐜1,2 = (±c, 0,0), 𝐜3,4 = (0, ±c, 0), 𝐜5,6 = (0, 0, ±c), 𝐜7−10 = (±c, ±c, 0),𝐜11−14 = (±c, 0, ±c), 𝐜15−18 =eq(0, ±c, ±c), Ql is the increment of gl attributable to the scalar source, and gl denotes the
local equilibrium distribution function. The second term on the right-hand side (r.h.s.) of
Eq. 14 represents the collision process, in which, after receiving the half-influence fromeq
the scalar source, the function gl gradually approaches gl at a time-rate defined by the collisional relaxation time 𝜏g (Guo et al. 2002). The equilibrium distribution function is givenby (Chopard et al. 2009)
]
[
)2 3
9(
eq
gl = El s 1 + 3𝐜l ∙ 𝐮 + 𝐜l ∙ 𝐮 − 𝐮 ∙ 𝐮 ,(15)
22
where s denotes the macroscopic scalar density and El denotes the weighing constants
E0 = 1∕3, E1−6 = 1∕18, and E7−18 = 1∕36. The source term in Eq. 14 is given by][
Ql = El qs + 3qs 𝐜l ∙ 𝐮 + 3rs 𝐜l ∙ 𝐅t ,(16)
where qs denotes the scalar source (or sink if negative), rs is the scalar mixing ratio definedas rs = s∕𝜌, and 𝐅t is the total force, including external forces and the pressure gradient
force. The quadratic terms in Eq. 15 and the third term in Eq. 16, which are not always considered in the literature, are included here to minimize spurious terms that would appear
in the macroscopic advection–diffusion equation reproduced from Eq. 14 (Chopard et al.
2009). Finally, using the zeroth-order moment of the distribution function, the scalar density is obtained as (Aursjø et al. 2017)

s=

18
∑

1
gl + qs .
2
l=0

(17)

By considering one-half of the source term here and in the collision term in Eq. 14, theoverall numerical accuracy is maintained at the second order both in time and space (Guoet al. 2002; Krüger et al. 2017).
The scalar diffusion resulting from unresolved eddies is incorporated via an SGS gradient–diffusion model as𝜏s𝛼 = −𝜅SGS

𝜕s
,
𝜕x𝛼

(18)

where 𝜏s𝛼 represents the SGS scalar flux and 𝜅SGS denotes the subgrid eddy diffusivity,
which is parametrized using the subgrid eddy viscosity and the subgrid turbulent Schmidt
number Sc as𝜅SGS = 𝜐SGS ∕Sc.

(19)

The effect of the SGS scalar flux is reflected in the simulations by including 𝜅SGS in thedefinition of the collisional relaxation time as13

T. Watanabe et al.

46

𝜅 + 𝜅SGS =

(
)
1
1
,
𝜏g −
3
2

(20)

where 𝜅 denotes the molecular diffusivity for the scalar. As a result, the present lattice
Boltzmann method for a scalar predicts the macroscopic scalar density that is governed by
the following filtered advection–diffusion equation
[
]
) 𝜕s
𝜕 (
𝜕s 𝜕u𝛼 s+=𝜅 + 𝜅SGS+ qs .
(21)
𝜕t
𝜕x𝛼
𝜕x𝛼
𝜕x𝛼

2.2.3 Simulation Set‑up
The numerical set-up for the flow simulation is the same as that used by Watanabe et al.
(2020). However, for the sake of completeness, it is described here along with the scalar
field. The simulation of a neutrally stratified turbulent flow within and above an idealized
uniform plant canopy is performed in a computational domain of 1024, 512, and 80 nodes
in the streamwise, lateral, and vertical directions, respectively. The grid spacing is 0.1 h
in all three coordinates, with h being the height of the plant canopy. Thus, the lowest 10
nodes represent the canopy layer, in which momentum and scalar sinks are distributed. The
airflow in the computationalis driven by a streamwise external force F ext , which()domain2−7
ext
is fixed at F ∕ 𝜌0 Δx∕(Δt) = 10 . The top boundary is a rigid free-slip wall at which
the scalar density is fixed to stop ∕s0 = 1, whereas the bottom boundary underlying the plant
canopy is a rigid wall, which exerts friction on the flow and absorbs the scalar through the
SGS fluxes of momentum and the scalar (defined later in Eq. 24), respectively. The horizontal boundaries are periodic for both the flow and scalar. The rigid top boundary for theflow and the rigid bottom boundaries of the flow and scalar fields are represented by meansof halfway specular reflection, in which distribution functions that stream into the domainfrom the wall are given by
(
)
(
)
∗
fijk x, y, zb , t + Δt = fijk
� x − icΔt, y − jcΔt, zb , t ,
(22a)

(
)
(
)
gl x, y, zb , t + Δt = g∗l� x − cl� x Δt, y − cl� y Δt, zb , t ,(22b)where zb is the height of the highest or lowest fluid node, and k and l denote inwards direc′′
tions opposite to the directions towards the wall (k and l ). This implies that the actual
wall boundaries are located one-half grid spacing beyond the highest or lowest fluid nodes(Lallemand and Luo 2003); all variables ( f , g, u, v, w, 𝜌, p, s) are collocated at each gridnode. The top boundary for the scalar distribution function is specified using the method of
Yoshino and Inamuro (2003) at the highest grid nodes (i.e. one-half grid spacing below the
top wall) to maintain the given constant scalar density at this height. The effects of plantelements on the flow and scalar fields are represented by an instantaneous drag force 𝐅dand a source of the scalar qs, parametrized as

13

𝐅d = −𝜌cd a|𝐮|𝐮,

(23a)

)
(
qs = −cs a|𝐮| s − sc ,

(23b)

47

Coherent Eddies Transporting Passive Scalars Through the Plant…

where cd and cs are, respectively, the drag coefficient and the scalar exchange coefficientof canopy elements, a is the leaf area density, and sc is the scalar density at leaf surfaces.
The canopy parameters are given as ah = 2 (i.e. the leaf area index LAI = 2), cd = 0.2,
cs ∕cd = 0.2, and sc = 0 to enable idealized simulations excluding any complexity arisen
from the canopy morphology and plant physiological factors. The frictional force and thescalar source from the underlying ground surface are given by

𝐅g = −𝜌b CMg ||𝐮b ||𝐮b ,

(24a)

)
(
qsg = −Csg ||𝐮b || sb − sg ,(24b)where CMg and Csg are the bulk transfer coefficients of momentum and the scalar, respectively, sg is the scalar density (set = 0) at the ground surface, and the subscript b indicates the values at the lowest nodes. The bulk transfer coefficients are evaluated assuming logarithmic profiles of the velocity and the scalar density with constant roughness
lengths for momentum ( z0g) and the scalar ( zsg), which are prescribed as z0g ∕h = 10−3 andzsg ∕z0g = 0.1, respectively, throughout the bottom boundary. The total force in Eq. 16 isthus given as 𝐅t = 𝐅ext + 𝐅d + 𝐅g − ∇p. The macroscopic velocity satisfying Eqs. 5 and23a is evaluated as (Guo and Zhao 2002)
∼

𝐮=
1
+
2

√

𝐮
∼
1
+ 12 cd a| 𝐮 |
4

,

(25)

∼

where 𝐮 denotes the macroscopic velocity calculated without the drag force from Eq. 5,
while the macroscopic scalar density satisfying Eqs. 17 and 23b is calculated froms=s̃
1 + 12 cs a|𝐮|

,

(26)

in which s̃ denotes the scalar density calculated from Eq. 17 without the scalar source(term. At the) lowest( nodes, )the expressions cd a and cs a in these equations are replaced by
cd a + CMg and cs a + Csg , respectively, to account for the friction and the scalar sourcefrom the ground. The subgrid turbulent Schmidt number is assigned a constant value of
Sc = 0.3, which is optimized to obtain scalar power spectra comparable with those fromthe Navier–Stokes model in the region well above the canopy.
The simulation is initialized using a uniform streamwise velocity component witha small amount of random perturbations in the vertical velocity component and a constant scalar density. The distribution functions are initially equilibrated with these macroscopic variables. After the turbulence statistics reach a steady state, three-dimensional
data pertaining to the macroscopic variables are collected at constant time intervals of
13.2h∕u∗; a total of 160 samples are used to calculate the statistics presented in the next
section. Time-sequential three-dimensional data are also collected every 1­ 03 timesteps(≈ 0.26h∕u∗) during the time periods between ± ­104 steps relative to the reference times,
which are set every ­105 steps in the first 3 × ­106 steps of the data-sampling period.
The mean values at each height are calculated as the ensemble mean of horizontally
averaged values obtained at a given height in each sample. In the following section,
an overbar (−) denotes the mean value, while a prime ( ′ ) indicates deviations from the
mean.

13

T. Watanabe et al.

48

3 Results and Discussion
The results are presented in dimensionless form using the friction velocity u∗, which isdefined at the canopy top by (Watanabe et al. 2020)()
u2∗ = F ext ztop − h ∕𝜌0 ,(27)
and the characteristic scale of the scalar s∗, estimated from()
s∗ = − s� w� + 𝜏s3 ∕u∗ ,a(28)

where the subscript a denotes a mean value averaged over the height range 2 ≤ z∕h ≤ 7.

3.1 Comparison with the Navier–Stokes Model
The performance of the lattice Boltzmann method is examined by comparing the resultswith those from a reference simulation by the established Navier–Stokes model. As thevalidation for the airflow field has been described elsewhere (Watanabe et al. 2020), thefocus here is the method’s ability to simulate a scalar field.

3.1.1 Mean Profiles
Figure 1 shows profiles of the mean scalar density, the resolved scalar variance, and thevertical scalar flux, calculated from the lattice Boltzmann and Navier–Stokes results. As aconstant scalar density is imposed at the top boundary, the mean density and the variancefrom both models increase sharply with height near the top boundary. This boundary condition also produces the vertically constant scalar flux seen above the canopy. There is anappreciable difference between the models in terms of the variance in the upper region of
the domain, which is probably attributable to the difference in the actual implementation of
the top boundary condition. In the lattice Boltzmann method, the downward-moving components of scalar distribution functions are specified at a half grid spacing below the topof the domain, whereas the downward SGS flux is evaluated exactly at the top boundary in
the Navier–Stokes model. As the leaf area is distributed uniformly within the canopy, andthe scalar density is also constant (= 0) at the surfaces of the canopy elements and ground,
there are no local maxima or minima in the distribution of the scalar sink/source; consequently, the mean scalar density, scalar variance, and scalar flux decrease monotonicallywith depth into the canopy. The counter-gradient fluxes are not simulated in the present
settings. Although these features (except for the constant flux) may not be typical for thescalar statistics of natural canopy turbulence, the figure indicates that the lattice Boltzmann
method reproduces the profiles simulated by the Navier–Stokes model, at least in the lower
half of the domain.

3.1.2 Budgets of the Scalar Variance and Flux
The budget equation for the resolved scalar variance is written as

13

Coherent Eddies Transporting Passive Scalars Through the Plant…

49

Fig. 1  Profiles of a mean scalar density, b resolved scalar variance, and c scalar flux calculated from thelattice Boltzmann method (LBM) (mark and line) and the Navier–Stokes (NS) model (line). Broken linesin (c) represent the SGS component of the scalar flux, and grey shading indicates the canopy layer in which
the leaf area is distributed homogeneously

(
𝜕
𝜕t

)
s�2
2

= 0 = −s� w�𝜕s 𝜕s�2 w� 𝜕s� 𝜏s3−−− cs a|𝐮|ss� − 𝜀s − Gdiv ,
𝜕z
𝜕z
𝜕z

(29)

where 𝜀s denotes the rate at which the resolved variance is converted to the SGS variance,
which is ultimately dissipated by the molecular diffusion, and Gdiv is a term that arises fromthe compressibility of the fluid, written as

13

T. Watanabe et al.

50

Gdiv =

1 𝜕ws�2 ss� 𝜕u𝛼
+
.
2 𝜕z
2 𝜕x𝛼

(30)

This term is included for completeness, although it is generally negligible in the present
results. The other four terms on the r.h.s. of Eq. 29 represent the rate of variance production from the interaction between turbulence and the mean scalar gradient, the resolved
turbulent transport, the SGS transport, and the loss of variance through the canopy sinks,
respectively. The( SGS dissipation
term for the Navier–Stokes model is calculated directly)��from 𝜀s = −𝜏 s𝛼 𝜕s ∕𝜕x𝛼 , whereas the sum of 𝜀s + Gdiv is evaluated as the residual of
Eq. 29 for the lattice Boltzmann model. The budget equation for the resolved scalar flux
can be written as𝜕s 𝜕s� w�2
𝜕s� w�
= 0 = −w�2 −
− Dsw
− Dsw
,
c
p
𝜕t
𝜕z
𝜕z

(31)

where the first two terms on the r.h.s. represent the flux production and the vertical transdenotes the decorrelation attributport by resolved turbulence, respectively. The term Dswcable to the canopy drag and sink��Dsw
c = cd a|𝐮|ws + cs a|𝐮|sw ,(32), which is evaluated as the residual of Eq. 31 for both models, represents
and the term Dswpthe other decorrelation effects caused by pressure and the SGS components.
Figure 2 shows profiles of each term in these budget equations calculated from theresults of the two models. It is clear that the lattice Boltzmann model reproduces the
production, transport, and dissipation or decorrelation of the scalar variance and flux
as well as the Navier–Stokes model. In Fig. 2a, both models predict that the scalar variance is locally in equilibrium for the region z∕h > 1.7, where the production is almost
balanced by the SGS dissipation. However, in the region below, the resolved and SGS
turbulent transport play a role in redistributing the variance from near the canopy top to
inside the canopy, where the scalar variance is dumped by the distributed scalar sinks.
Figure 2b indicates that the scalar flux is produced by the interaction between the variance of the vertical velocity component and the mean scalar gradient, and then lost by
the decorrelation processes caused by the canopy sink and the total effect of the pressure and SGS components. The turbulent transport plays a more important role than thevariance budget in redistributing the flux from above the canopy to the interior, resulting in a slower equilibration (attained at z∕h ≳ 4 ) than the scalar variance. The production terms in both budgets peak sharply at the canopy top, but decrease quickly withdepth into the canopy as the vertical scalar flux and the variance of the vertical velocity
component diminish. Close to the canopy floor, the largest sources of the scalar varianceand flux are provided by the turbulent transport, represented by the third-order mixed
moments, demonstrating the importance of the non-local scalar transport in this region.
The significant contribution of the canopy terms to both budgets is the result of theparametrizations of the scalar sink and drag force by the canopy elements, whereby they
are dependent on the instantaneous scalar density and velocities (Eq. 23). If a constantsink/source strength and a constant drag force were adopted to represent the plant canopy, the canopy terms would vanish from both budgets. Hence, the fluctuations of sink/
source strength should not be overlooked when considering the budgets of the secondorder moments of scalar turbulence in plant canopies.

13

Coherent Eddies Transporting Passive Scalars Through the Plant…

51

Fig. 2  Profiles of the terms in the budget equations of a scalar variance and b scalar flux simulated by thelattice Boltzmann (lines) and Navier–Stokes (marks) models. Ps: production, Tt: turbulent transport, Dc:
canopy sink, Dp: decorrelation by pressure and SGS components, Tsgs: SGS transport, Dsgs: loss by SGS
production. The Navier–Stokes results shown here have been smoothed with a 1–2–1 filter

3.1.3 Power Spectra of the Scalar
Figure 3 shows horizontal one-dimensional power spectra of the scalar fluctuations atselected levels. The simulated spectra generally follow the −5∕3 power law at levels
above the canopy but have slightly steeper gradients within the canopy, consistent with

13

52

T. Watanabe et al.

Fig. 3  One-dimensional power spectra of the scalar variance: a streamwise and b lateral spectra calculated
for four different heights of z∕h = 3.05, 2.05, 1.05, and 0.55. Lines represent the lattice Boltzmann results,
and marks denote the results of the Navier–Stokes model; Lx and Ly are the streamwise and lateral dimensions of the computational domain; 𝜆x and 𝜆y are the streamwise and lateral wavelengths. The topmost spectra are aligned to the vertical axis, and others are displaced vertically for clarity (multiplied by 0.5, 0.25,
and 0.125, respectively)

previous observations and experiments (Maitani and Seo 1985; Poggi et al. 2011; Katul
et al. 2013). The spectra from the two models overlap with each other for most wavenumbers. The Navier–Stokes spectra show a signature of the numerical aliasing at thelargest wavenumbers. Although the horizontal derivatives are calculated by the pseudospectral method with the two-thirds dealiasing rule, the Navier–Stokes model used here
may be subject to dispersion errors from the vertical finite-differencing, especially inthe regions of sharp scalar fronts caused by intense vertical motions. Adopting a lessdispersive advection scheme and/or a more effective SGS diffusion scheme that selectively damps such errors may therefore be desirable. The lattice Boltzmann model, however, tends to attenuate the scalar fluctuations for wavelengths shorter than about 5Δ;
this considerably enhances the numerical stability, but also requires a finer grid systemthan the Navier–Stokes model to attain a given resolution. In the small wavenumber
range, the lateral spectra above the canopy have peaks around a wavenumber of 0.06h∕𝜆y , whereas the streamwise spectra gradually increase with decreasing wavenumber,indicating streamwise structures in the periodic computational domain.

3.1.4 Two‑point Autocorrelation Functions of the Scalar Density
The spatial structure of the coherent eddies that induce energetic perturbations at agiven reference height zr can be captured by the two-point correlation function between
the flow or scalar variables obtained at an arbitrary height and at the reference height.
This function can be written as

13

53

Coherent Eddies Transporting Passive Scalars Through the Plant…

(

Rab rx , ry , z, zr

)

)
(
) (
a� x + rx , y + ry , z b� x, y, zr=[( )]1∕2a�2 (z) b�2 zr(33)

where a and b denote the variables, and rx and ry denote the streamwise and lateral separation, respectively. Figure 4 shows the spatial distribution of the scalar’s autocorrelation
function Rss for a reference height just above the canopy at zr ∕h = 1.05. The figure reveals
that the scalar perturbations at the canopy top are typically associated with spatial distributions of the scalar density that are streamwise-elongated (Fig. 4a), relatively narrowin the lateral direction (Fig. 4b), vertically tall, and tilted downstream (Fig. 4c), similar
to a previous numerical study by Su et al. (2000) and a field observation by Dupont and
Patton (2012). Although the highest correlations are confined within a narrow, inclinedcore region residing near the canopy top, a peripheral positive correlation with Rss > 0.1
extends beyond the figure both downstream and upstream. Downstream of the core region,
the positive correlation spreads towards far above the canopy, whereas a long tail of positive correlation continues within the canopy upstream of the core region. The streamwiseelongated and downstream-tilted features of Rss are similar to those of Ruu reported by
Watanabe et al. (2020). However, the enhanced vertical gradient of the mean scalar density(Fig. 1a), which is caused by the imposed constant top boundary, may somewhat exaggerate the vertical extent of the scalar distribution.
The above comparisons confirm that the lattice Boltzmann method is as capable as theNavier–Stokes model of simulating the scalar turbulence within and above a plant canopy.
Although direct comparisons with observations from the actual plant canopies are difficult
as an idealized uniform canopy is considered here, the models reproduce well-known general characteristics of scalar turbulence described in the Introduction. In the following, thecharacteristics and spatial structures of the turbulent eddies that cause energetic perturbations in the scalar density and promote the vertical scalar transport through the canopy are
investigated using the lattice Boltzmann results.

3.2 Characteristics of Eddies Associated with Scalar Variance and Flux Near
the Canopy Top
3.2.1 Two‑point Cross‑correlation Functions Between Velocity and Scalar
The two-point cross-correlation functions between the velocity components (u and w ) at
arbitrary heights and the scalar density ( s ) at the reference height are shown in Fig. 5.
Both velocity components have high correlations with the scalar density, with the maximum magnitude of the correlation function being greater than 0.5. The spatial distribution of Rus exhibits a very large, streamwise-tilted pattern, and that of Rws is verticallycoherent and horizontally compact. These patterns of Rus and Rws closely resemble thespatial structures of Ruu and Rww reported by Watanabe et al. (2020) and indicate thatthe eddies inducing the predominant perturbations in the streamwise and vertical velocity components are also responsible for the predominant scalar perturbations near the
canopy top. Because sweep motions generally dominate in the TKE near the canopy top
and in the scalar flux under the current settings (see the quadrant analysis shown below),′
the most significant scalar variations should be positive perturbations ( s > 0 ) caused′′by eddies with u > 0 and w < 0 (Fig. 5c). As demonstrated by Watanabe (2004), whena high-speed downward (sweep) motion carrying air of high scalar concentration from

13

54

T. Watanabe et al.

Fig. 4  Spatial distribution of two-point autocorrelation functions of the scalar density, Rss, in a horizontal
cross section at z∕h = 1.05, b vertical–lateral cross section, and c vertical–streamwise cross section. Line
contours with intervals of 0.1 denote the lattice Boltzmann results, coloured shading represents the Navier–Stokes results

above affects the canopy top, the scalar density around the impact position increases
sharply, producing a ramp signal in the scalar time traces. As the sweep motion translates downstream faster than surrounding air masses with lower scalar density, a downstream-tilted scalar microfront is generated in front of the sweep motion if properly

13

Coherent Eddies Transporting Passive Scalars Through the Plant…

55

Fig. 5  Vertical–streamwise distribution of two-point cross-correlation functions of a streamwise velocity
component ( Rus) and b vertical velocity component ( Rws), referenced to the scalar density at a height justabove the canopy ( z∕h = 1.05). In c, these functions are also displayed as vectors overlaid on the map of thetwo-point autocorrelation function of scalar density shown in Fig. 4c. The interval of line contours is 0.1

aligned with ejection motions (Gao et al. 1989; Fitzmaurise et al. 2004; Finnigan et al.
2009). After the sweep motion penetrates the canopy layer, the translation speed isdecelerated by the drag of canopy elements, and a weak reverse flow is induced near
the ground (Fig. 5a), as pointed out by Su et al. (2000) and Shaw et al. (2013). This
may enhance the long tail of the scalar trailing upstream of the impact position. Thus, acombination of the sheared mean velocity above the canopy and the retarded translation
speed inside the canopy causes the overall structure of the scalar to tilt downstream. The

13

T. Watanabe et al.

56

long tail of the same structure may cause a scalar-sweep event near the canopy floor, as
will be discussed in Sect. 3.3.

3.2.2 Quadrant Analysis
The joint probability density function (p.d.f.) of s and w for four different heights near the
canopy top is displayed in Fig. 6. The figure also shows bold-line contours of the third-ordertruncated Gram–Charlier distribution (Antonia and Atkinson 1973; Nakagawa and Nezu
1977), given by[( )]∑m 𝜕 p+q G ̂s, w
̂
( )
( )
p+q pqp ̂s, w
̂ = G ̂s, w
̂ 1+
,
(−1)
(34)
p
q
p!q!
𝜕̂s w
̂
p+q=3
′

′

̂ = w ∕𝜎w are the variables normalized by the corresponding standard
where ̂s = s ∕𝜎s and w
̂ q ) ( p, q ∈ {0,1, 2,3}; p + q = 3) are the normalized third-orderdeviations and (mpq = ̂sp w
��

Fig. 6  Joint p.d.f.’s of normalized fluctuations of the scalar density and the vertical velocity component atfour different heights: a z∕h = 0.55, b 1.05, c 1.25, and d 1.55. The coloured shading and dashed line contours represent the simulation results, while the bold line contours denote the third-order truncated Gram–
Charlier distribution

13

57

Coherent Eddies Transporting Passive Scalars Through the Plant…

central( ) moments, or equivalently the normalized third-order cumulants. The function
G ̂s, w
̂ represents the two-dimensional Gaussian
]
[
( )
̂+w
̂2
̂s2 − 2Rsŵsw
1
,
exp −
G ̂s, w
̂ =
(
)
√
(35)
2 1 − R2sw
2𝜋 1 − R2sw

̂ denoting the correlation coefficient between s and w . As the covariancewith Rsw = ̂sw
s′ w′ is negative in the present simulation, each quadrant of the coordinate is classified asfollows:
� > 0) Outward interaction,
Quadrant I ( �s > 0, w
� > 0) Ejection,
Quadrant II ( �s < 0, w
� < 0) Inward interaction,
Quadrant III ( �s < 0, w
� < 0) Sweep.
Quadrant IV ( �s > 0, w′′

Figure 6 shows that the joint p.d.f. at each height has a peak in quadrant II and a gently
sloping skirt extending towards quadrant IV, consistent with previous analyses (Coppin
et al. 1986; Maitani and Shaw1990; Dupont and Patton 2012) and confirming the wellknown fact that sweeps dominate over ejections in the vertical turbulent transport at theseheights near the canopy top (Finnigan 2000). The truncated Gram–Charlier distribution
approximates the joint p.d.f. sufficiently far above the canopy ( z∕h = 1.55). However,near and inside the canopy layer, the approximation becomes slightly worse, indicating
the influence of higher-order deviations from the Gaussian. At these heights ( z∕h = 1.25,
1.05, and 0.75), the actual p.d.f. within a range of small positive values of the ordinate
� < 1) spreads in the negative region of the abscissa towards the minimum of the(0 < wscalar perturbation at each level. A similar tendency has also been observed for the joint′
p.d.f. of temperature and w measured in a deciduous forest (Maitani and Shaw 1990), as′′
well as that of u and w measured within the canopy layer of a deciduous forest (Shaw et al.
1989; Maitani and Shaw 1990) and a model canopy in a wind tunnel (Raupach 1981). Such
non-Gaussian behaviour near the canopy top is likely to be related to the non-local transport caused by the ejection motions that pump up the air, which has depleted (or enriched)scalar density and reduced streamwise momentum, from near the bottom to the top of the
canopy. This point will be discussed in more detail later in the next section with respect tothe conditional sampling analysis for scalar ejection events near the canopy floor.
Admitting the third-order Gram–Charlier distribution as a reasonable approximation, it′
is worth examining the third-order moments of s′ and w in relation to the eddy motions
in each quadrant (Raupach 1981; Katul et al. 1997a, b, 2018; Cava et al. 2006). From thisidealized simulation, pure information can be obtained on how the non-local transports are
related with the second moments without any complexity arisen from canopy morphology, biological factors, and the atmospheric conditions, which are unavoidable in natural
canopies. Using the third-order Gram–Charlier distribution function, the differences inthe second moments calculated for each quadrant are analytically related to the third-ordermoments by (Nakagawa and Nezu 1977; Raupach 1981)� � � �
�
�
��
1 � �
̂sw
̂ I − ̂sw
̂ III = √3 m21 + m12 − Rsw m30 + m03 ,
6 2𝜋

(36a)

13

T. Watanabe et al.

58

� �
� �
�
�
��
1 � �
̂sw
̂ IV − ̂sŵ II = √3 m21 − m12 − Rsw m30 − m03 ,
6 2𝜋

(36b)

� � � �� �� �2̂ŝs I − ̂ŝs II − ̂ŝs III + ̂ŝs IV = √ m30 ,
3 2𝜋

(37a)

� � � �
� �
� �
2
w
̂w
̂ I+ w
̂w
̂ II − w
̂w
̂ III − w
̂w
̂ IV = √ m03 .
3 2𝜋

(37b)

The third-order mixed moments, m12 and m21, which are related to the vertical transport ofthe variance (s�2 ) and the covariance s′ w′ , respectively, are thus written in terms of the second moments in each quadrant as√ �� � � �
� �
� � �
̂ I − ̂sw
̂ II − ̂sŵ III + ̂sw
̂ IV +
m21 = 2𝜋 ̂sw√�� � � �
� �
� � �
2𝜋
Rsw ̂ŝs I − ̂ŝs II − ̂ŝs III + ̂ŝs IV ,
2

(38a)

√
�� � � �√ �� � � �
� �
� � �
� �
� � �
2𝜋
̂ I + ̂sw
̂ II − ̂sŵ III − ̂sŵ IV +
̂w
̂ I+ w
̂w
̂ II − w
̂w
̂ III − w
̂w
̂ IV .
Rsw w
m12 = 2𝜋 ̂sw
2(38b)Figure 7a shows that the third-order moments estimated by these equations closely
match the actual moments calculated directly from the simulated turbulence field.

Fig. 7  Vertical profiles of a the normalized third moments of scalar density and vertical velocity components and b the ratio of m12 to m21. Marks denote the values calculated directly from the lattice Boltzmann
results, and lines represent the estimations from Eqs. 37 and 38

13

Coherent Eddies Transporting Passive Scalars Through the Plant…

59

Equation 37 implies that the skewness can be viewed as the difference between the variance of positive and negative perturbations. Similarly, noting that the first terms on the
r.h.s. of Eq. 38 are generally more important than the second terms (representing the effectof skewness), the moments m12 (or m21) essentially represent the difference between the′′contributions of eddy motions with positive and negative perturbations of w (or s ) to thecovariance. Thus, non-Gaussian motions such as coherent eddies developing in the canopy turbulence, which cause asymmetry in the variances and covariances associated withdownward and upward motions (especially sweeps and ejections), is responsible for thelarge contributions of the turbulent transport terms to the budgets of s′2 and s′ w′ (Fig. 2).
These descriptions are consistent with previous studies, demonstrating that the third-order
Gram–Charlier distribution predicts many aspects of the joint p.d.f., including the imbalance between the contributions of sweeps and ejections to the transport of momentum andscalars and its link to the third-order moments measured in the surface layers, including
the roughness and canopy layers (Raupach 1981; Shaw et al. 1983; Katul et al. 1997a, b,2018). Based on field observations, previous studies have also shown that the( moments)( m)21
̂
s
w
̂
−
̂sw
̂ II
imbalance
betweensweepsandejectionsand m12 are directly(related to theIV)when the term Rsw m30 − m03 is neglected and a linear relationship between m21 and m12
is assumed in Eq. 36b (Cava et al. 2006; Katul et al. 2018). However, such simplifications
are not adequate for (the present)simulation results, especially above the canopy. This isbecause the term Rsw m30 − m03 is not negligible—m30 has a large value, which may be
attributable to the constant scalar density imposed at the top boundary—and the contributions to the covariance from quadrants I and III are not linearly related with those from
quadrants II and IV thereby causing the ratio m12 ∕m21 to vary with height (Fig. 7b).

3.3 Coherent Eddy Transporting Scalars to/from Deep Inside the Canopy
As indicated by Fig. 2, the major source of the variance and flux of the scalar near the canopy floor is the non-local turbulent transport driven by large-scale coherent eddies, whichcause a higher-order non-Gaussian signature in the joint p.d.f. near the canopy top (Fig. 6).
In this subsection, the structure of canopy-penetrating eddies that drive the non-local transport of scalars between near the canopy floor and above the canopy is examined.

3.3.1 Instantaneous Snapshot
Figure 8 shows snapshots of the horizontal distributions of velocity fluctuations and thescalar sweeps and ejections at a height near the canopy floor ( z∕h = 0.15). In contrast to
the streaky patterns that dominate the velocity distributions above the canopy (Watanabe
2004, 2009), cellular patterns are noticeable in the distribution of the vertical velocity component (Fig. 8a), of which positive and negative areas correspond to the zones of convergence and divergence, respectively, of horizontal velocity components. As indicated by the
TKE budgets reported in previous numerical studies (Dwyer et al. 1997; Watanabe et al.
2020), the primary source of TKE near the canopy floor is the pressure transport, that is,
the kinetic energy of vertical motions caused by the work of pressure. However, as verticalmotions are blocked by the solid ground surface, horizontal divergence and convergence
are induced (corresponding to the downdrafts and updrafts, respectively), thus producing
the cellular patterns. This view is supported by the fact that, in the absence of strong shear
deep inside the canopy, the correlation between the pressure and the normal strain rate,
which redistributes energy among velocity components, provides the major source for the

13

60

T. Watanabe et al.

Fig. 8  Instantaneous distributions in a horizontal plane at z∕h = 0.15 of a velocity fluctuations and b scalar flux in sweep and ejection events. Vectors represent fluctuations of the horizontal velocity components�(u� ∕u∗, v� ∕u∗(), and) colour shading indicates the vertical velocity component w ∕u∗ in (a) and the scalar
fluxes s� w� ∕ u∗ s∗ by sweep (blue) and ejection (red) motions in (b) with the sign of the ejection flux
reversedvariance of horizontal velocity components at the expense of the vertical velocity component near the canopy floor (figure not shown). The tight connection between the pressure
fluctuations and horizontal velocity components near the canopy floor has been identified
in field measurements (Shaw and Zhang 1992; Zhang and Amiro 1994) and in numericalsimulation results (Su et al. 2000; Shaw et al. 2013). In this velocity field, the scalar sweep
and ejection events are distributed according to the vertical velocity component (Fig. 8b);the sweep events can be seen in regions of strong downdrafts and the ejection events occur
in the cellular convergence zone surrounding the downdrafts.

13

Coherent Eddies Transporting Passive Scalars Through the Plant…

61

3.3.2 Ensemble‑averaged Structure of the Scalar Sweeps and Ejections Near
the Canopy Floor
The mean spatial structure of such canopy-penetrating eddies can be deduced by compositing the conditionally sampled variables as�∑n ∑mj � � ic xr + rx , yir + ry , z, tj + 𝜏
�
�
j=1
i=1
�
⟨c ⟩ rx , ry , z, zr , 𝜏 =
,
(39)
∑n
j=1 mjwhere ⟨c′ ⟩ denotes an ensemble average of any variable c′ sampled at time tj + 𝜏 (with 𝜏
representing a time lag) around detection positions (xri , yir , zr ) (i = 1, … , mj ), where thedetection conditions are satisfied in the j-th sample obtained at time tj , rx and ry are thestreamwise and lateral displacements from the detection position, respectively, n is thenumber of samples, and mj is the number of detected events in the j-th sample. The heightof the detection position is set close to the canopy floor at zr ∕h = 0.15, and the conditions
for the sweep and ejection events are set( as follows:)�
Scalar sweep: w′ < 0 and s w� ∕( 𝜎s 𝜎w) < −2,�Scalar ejection: w′ > 0 and s w� ∕ 𝜎s 𝜎w < −2.
Figures 9 and 10 show vertical slices of the no-lag (𝜏 = 0) composite results for sweep
and ejection events. Figure 9 indicates that the scalar sweep events near the canopy floor
most often occur near the trailing (upstream) ends of very large sweep motions passing
over the canopy. Both the streamwise velocity component and the vertical velocity component exhibit a very large pattern spreading above the canopy and tilting downstream. A
weak negative pressure perturbation is centred inside the canopy, between the downdraft ofthe main sweep motion and a weak updraft slightly upstream from the detection position.
The pattern of scalar density exhibits a downstream-tilted, tall structure above the canopy
and a long tail inside the canopy, which is similar to the distribution of the two-point auto′
correlation of s shown in Figs. 4 and 5. This similarity can be interpreted as the large-scale
sweep motions, which cause the predominant perturbations in the scalar density near the
canopy top, also inducing scalar sweep events near the canopy floor at their trailing ends.
The scalar ejection events are associated with more complex patterns of the streamwise
velocity component (Fig. 10a). Near the canopy top, above the detection position, there isa negative perturbation of u, which continues inside the canopy and shifts downstream as
it approaches the canopy floor. A shallow positive perturbation is located on the canopy
floor just upstream of the detection position, and another one appears around rx ∕h = 3 − 7,
above which a downstream-tilted positive perturbation extends from near the canopy top
to approximately z∕h = 5. Further downstream, a negative perturbation of the streamwise
velocity component extends above the canopy. The patterns in the vertical velocity component (Fig. 10b) are characterized by an updraft/downdraft pair in which the updraft has
the greater magnitude. The shape of the updraft/downdraft pair is narrow and verticallycoherent, in contrast to the scalar sweep, which has a larger structure tilting downstream.
The maximum speed of the updraft occurs exactly at the height of the canopy top above thedetection position, while that of the downdraft appears downstream and is slightly liftedfrom the canopy top. In Fig. 10c, a negative pressure core, which is more intense than
the sweep case, is also centred at the canopy top between the paired updraft and downdraft. Figure 10d illustrates that, near the detection position, the coherent updraft carries
air with lower scalar density from near the canopy floor to above the canopy. On the downstream side, a downstream-tilted area with higher scalar density, similar to that seen in the

13

62

T. Watanabe et al.

Fig. 9  Vertical–streamwise cross-sections of the ensemble-averaged structure of eddies causing scalar�sweep events at a height close to the canopy floor ( zr ∕h �= 0.15
� ): a streamwise velocity� component ⟨u ⟩∕u∗,b vertical velocity component ⟨w� ⟩∕u∗, c pressure ⟨p� ⟩∕ 𝜌u2∗ , and d scalar density ⟨s ⟩∕s∗ with vectors of(⟨u� ⟩∕u∗, ⟨w� ⟩∕u∗). Contour intervals are a 0.2, b 0.1, c 0.2, and d 0.3

streamwise velocity component, lies above an area of relatively lower scalar density that isaccompanied by weak positive perturbations of the vertical velocity component.
Horizontal slices of the composite velocity fluctuations at the canopy top ( z∕h = 1.05)
also show that the downdraft region in the sweep case shifts downstream from the13

Coherent Eddies Transporting Passive Scalars Through the Plant…

63

Fig. 10  Same as Fig. 9, but for scalar ejection events

detection position, whereas the updraft centre in the ejection case is located exactly overthe detection position (Fig. 11). The most active parts of the eddy structures are streamwise-elongated (for sweeps) or horizontally compact (for ejections), rather than havingtransversally coherent forms. Downstream of both structures, there are streamwise-elongated regions of high- or low-speed streamwise velocity components, which may either
be caused by the eddies or provide a place for the eddies to develop.

13

64

T. Watanabe et al.

Fig. 11  Horizontal cross-sections at a height just above the canopy (z∕h = 1.05) of the ensemble-averaged
structures of eddies causing a scalar sweep and b scalar ejection events at a height close to the canopy floor( zr ∕h = 0.15). Vectors represent the horizontal velocity components (⟨u� ⟩∕u∗, ⟨v� ⟩∕u∗), and coloured shading indicates the vertical velocity component ⟨w� ⟩∕u∗3.3.3 Temporal Development of Canopy‑penetrating Eddies
Time sequences of the time-lagged composite of velocities and scalar density are presented
in Fig. 12 (for the scalar sweeps) and 13 (for the scalar ejections), in which only the lower
half-domain with a narrower horizontal range is shown. Figure 12 clearly demonstrates
that the scalar sweep events detected near the canopy floor originate from large-scale
sweep motions coming from aloft and penetrating the canopy. At the earliest time in the
figure, the centre of the positive scalar perturbation occurs at a height of around 2h. Conveyed by the sweeping motion, the scalar’s positive perturbation approaches the canopy
and is stretched diagonally to the downstream-upward direction. After the strongest part ofthe sweeping motion has passed over the detection position, the highest part of the scalar
perturbation touches down on the canopy floor, causing a scalar sweep event to occur at
this moment. The pattern of scalar perturbation above the canopy then becomes further
stretched downstream, with a long tail left behind inside the canopy for a similar reason as
discussed for Fig. 5.
Figure 13 shows that the scalar ejection events near the canopy floor are also triggered by sweep motions from aloft, although the process is not as straightforward as13

Coherent Eddies Transporting Passive Scalars Through the Plant…

65

Fig. 12  Vertical–streamwise cross-sections of time-lagged composite structures for eddies causing scalarsweep events at a height close to the canopy floor ( zr ∕h = 0.15). Time lags are a − 2.34, b − 1.56, c – 0.78,
d 0, and e + 0.78 in units of h∕u∗ . Coloured shading and line contours of an interval of 0.3 indicate the scalar density ⟨s� ⟩∕s∗ and vectors represent (⟨u� ⟩∕u∗, ⟨w� ⟩∕u∗)for the sweep case. The earliest time in the figure shows a build-up of air with low
scalar density, which is associated with low-speed ejective motions with large spatial
scales (Fig. 13a). The lowest scalar density is seen near the canopy floor upstream
from the detection position. While a massive sweep motion carrying high scalar densities from aloft approaches over the low-speed air of low scalar density (Fig. 13b), an
updraft is generated upstream of the motion, and a reverse flow strengthens inside thecanopy beneath the motion (caused by the pressure, as discussed in Sect. 3.3.1). As the
direction of velocity perturbations prior to the detection time is inclined so as to draw′
energy from the mean shear (i.e. ⟨u ⟩⟨w′ ⟩ < 0 ; Fig. 13c), the downdraft of the sweep

13

66

T. Watanabe et al.

Fig. 13  Same as Fig. 12, but for scalar ejection events

motion itself and the induced updraft are both intensified until they reach a maximum
shortly after the sweep motion passes over the detection position. The resulting coherent
updraft pumps up air with the lowest scalar density to above the canopy, causing large
negative fluctuations of the scalar density near the canopy top (Fig. 13d), as mentioned
earlier regarding the non-Gaussian nature of the joint p.d.f. (Sect. 3.2.2). Unlike the scalar sweeps, there is only a marginal delay between the occurrence of maximum updrafts
near the canopy top and the scalar ejection events near the canopy floor. The pattern ofejected lower scalar density is then lifted further by the remaining updraft and stretchedby the mean shear (Fig. 13e).
As the conditional sampling is triggered by the scalar-transporting events at a level
close to the canopy floor, the largest perturbation in the composited scalar density naturally appears near the canopy floor for both the sweep and ejection events (Figs. 9d and10d). The maximum values of the downward velocity component and pressure are also13

Coherent Eddies Transporting Passive Scalars Through the Plant…

67

seen inside the canopy for the sweep case (Fig. 9b, c). However, the largest magnitudes of
the upward velocity component and the negative pressure perturbations in the ejection caseare observed at the level of the canopy height (Fig. 10a, b). The vertical coherency of thesevariables is more pronounced than for the sweep case, in which the pattern of the verticalvelocity component is tilted downstream. The intensity of pressure signals is much larger
for the ejection case. These facts imply that the scalar ejection events are caused aroundthe time of maximum development of the coherent eddy generated near the canopy top(Fig. 13), and that pressure plays a more important role in inducing the vertical motions
near the canopy floor in the ejection case.
Finally, it is worth mentioning that the eddy structure obtained by compositing the scalar ejection events near the canopy floor (Fig. 10) resembles that obtained by compositing the prominent updraft/downdraft pairs detected at the canopy top by Watanabe (2009).
This previous study also identified the complex pattern of the streamwise velocity component, the updraft/downdraft pair located upstream and downstream of a vertically coherent
negative pressure core, respectively, and the downstream-tilted positive and negative perturbations of the scalar density. Watanabe (2009) attributed the structure to an ensembleaveraging of vortical eddies that developed near the canopy top. Although some differences are appreciable, such as the magnitude of the downdraft being greater than that of
the updraft and the positive pressure perturbation being located downstream of the negative pressure core (which are attributable to the completely different detection schemes
and positions adopted), the similarity of the composite structures is quite remarkable. Thissimilarity indicates that the non-local scalar ejections from the canopy floor are driven by
those coherent eddies that cause predominant perturbations of the vertical velocity component near the canopy top.
As for the mechanism of the coherent eddies characterizing the turbulence near the
canopy top, Finnigan et al. (2009) explained them as pairs of counter-rotating vortices
developed through the deflection of initially two-dimensional Stuart vortices originatedfrom the inflection instability of the mean wind profile. This formation process was confirmed by Gavrilov et al. (2013). In contrast, Bailey and Stoll (2016) suggested that thetwo-dimensional roller-like structures themselves remain to contribute the Reynolds stressnear the canopy even in a fully developed turbulence. A reasonable next step would thus be
to explore whether these eddies directly induce the non-local scalar transport to/from thecanopy floor or if additional mechanisms are involved.

4 Conclusions
We have described a double-distribution-function lattice Boltzmann model for simulating apassive scalar field in a neutrally stratified turbulent flow within and above a plant canopy
using the LES approach. The model is composed of a central-moment-based multi-relaxation-time model for the airflow and a single-relaxation-time model for a scalar. In the simulation of scalar turbulence within and above a homogeneous plant canopy, the performanceof the model was shown to be comparable to that of a conventional LES model based on
the Navier–Stokes equations and a scalar advection–diffusion equation. Therefore, users
can choose either models depending on the available computer resources to perform efficient simulations of the dispersion of passive scalars within and above plant canopies. If
its larger memory consumption is acceptable, the lattice Boltzmann method described here

13

68

T. Watanabe et al.

would be a suitable choice for graphics-processing-unit-based computers because of its
local and continuous memory access.
In the present simulation settings, in which a passive scalar is transferred downward fromthe atmosphere (i.e. a top-down scalar) by a neutrally stratified turbulent flow, the large fluctuations in scalar density that occur near the canopy top are predominantly caused by sweep
motions that have a spatial scale much larger than the canopy height. The downstream-tilted
spatial structures of the scalar density continue from inside the canopy up to at least a few
times the canopy height as a result of vertical differences in the translation velocity of thesweep motions in conjunction with the sheared mean wind speed. The sweep motions further
penetrate the canopy and make a significant contribution to the scalar flux within the canopy.
These simulated characteristics of the scalar turbulence are consistent with observations fromwithin and above plant canopies.
The scalar sweep events near the canopy floor are a direct consequence of the penetrationof such large-scale sweep motions, which push down the tailing part of the downstream-tilted
structure of the scalar density. By contrast, the scalar ejection events near the canopy floor
are driven by updrafts of coherent eddies generated near the canopy top under the conditions
that sweep motions coming from aloft approach existing wide regions of low-speed ejective
motions. The updrafts are intensified and tilted by the mean shear until they are almost inphase at levels from near the canopy floor to approximately twice the canopy height at theirmaximum strength, when scalar ejection events are induced near the canopy floor.
These coherent eddies contribute to the non-local transport of the scalars emitted/absorbedby the underlying soil surface to/from above the canopy. A signature of such a direct transport through the canopy is detectable even above the canopy as a skewed deviation of thejoint p.d.f. of the scalar and vertical velocity fluctuations from the third-order Gram–Charlier distribution. The sweeping downdrafts generally make a larger contribution to the flux
than the ejecting updrafts in the present simulation. However, if the turbulence velocity field isunchanged, the updrafts generated near or above the canopy top play a more significant role intransporting those scalars that have their primary sources at the soil surface (e.g., respiratory­CO2, ­CH4, volatile organic matter). The resulting change in the contributions of downdrafts
and updrafts to the vertical scalar flux also impacts the direction and magnitude of the turbulent transport of the scalar variance and flux, as indicated by an analysis using a third-order
Gram–Charlier approximation of the joint p.d.f. This analysis reinforces the importance oftaking the spatial distributions of the scalar source/sink within the canopy, as well as the structures and causal mechanism of coherent eddies, into consideration when modelling the nonlocal transport of scalars.
The characteristics of the coherent eddies described here were deduced from a single caseof a neutrally stratified flow in a relatively sparse plant canopy. The influences of the canopy
morphology and the thermal stratification on their structures and resulting non-local transports are left as future issues. The lattice Boltzmann model presented herein can be applied tothermally non-neutral cases with marginal modifications to the code and would be feasible forfurther investigations on the non-local transport phenomena in wider conditions.
Acknowledgements The authors would like to thank Gabriel Katul (Duke University) for his insightful suggestions on the scalar turbulence and Yuta Hasegawa (Japan Atomic Energy Agency) for invaluable discussions on the computational aspects of the lattice Boltzmann method. This study was funded, in part, by JSPS
KAKENHI (Grant Number JP20K04057) and by the Grant for Joint Research Program of the Institute of
Low Temperature Science (ILTS), Hokkaido University (Grant Numbers 19G007, 20G008). The reference
simulation was performed on the supercomputing system (NEC SX-9 as of 2012) at the Computer Centre
for Agriculture, Forestry and Fisheries Research Information Technology Centre, MAFF, Japan, during a

13

Coherent Eddies Transporting Passive Scalars Through the Plant…

69

joint research programme between the Forestry and Forest Products Research Institute (PI: Yasuko Mizoguchi) and ILTS. Thanks are also due to Stuart Jenkinson (Edanz Group) for English language editing.
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 longas 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 the article’s Creative Commons licence, unless indicated otherwise in a credit line to thematerial. If material is not included in the article’s Creative Commons licence and your intended use is not
permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directlyfrom the copyright holder. To view a copy of this licence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/.

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
Antonia RA, Atkinson JD (1973) High-order moments of Reynolds shear stress fluctuations in a turbulent boundary layer. J Fluid Mech 58:581–593
Aursjø O, Jettestuen E, Vinningland JL, Hiorth A (2017) An improved lattice Boltzmann method for
simulating advective–diffusive processes in fluids. J Comput Phys 332:363–375
Bergstöm H, Högström U (1989) Turbulent exchange above a pine forest. II organized structures.
Boundary-Layer Meteorol 49:231–263
Bailey BN, Stoll R (2016) The creation and evolution of coherent structures in plant canopy flows andtheir role in turbulent transport. J Fluid Mech 789:425–460
Brunet Y (2020) Turbulent flow in plant canopies: historical perspective and overview. Boundary-Layer
Meteorol 177:315–364
Cava D, Katul GG, Scrimieri A, Poggi D, Cescatti A, Giostra U (2006) Buoyancy and the sensible heat
flux budget within dense canopies. Boundary-Layer Meteorol 118:217–240
Chen S, Doolen GD (1998) Lattice Boltzmann method for fluid flows. Annu Rev Fluid Mech 30:329–364
Chopard B, Falcone JL, Latt J (2009) The lattice Boltzmann advection–diffusion model revisited. Eur
Phys J Spec Topics 171:245–249
Coppin PA, Raupach MR, Legg BJ (1986) Experiments on scalar dispersion within a model plant canopy. Part II: An elevated plane source. Boundary-Layer Meteorol 35:167–191
Deardorff JW (1971) On the magnitude of the subgrid scale eddy coefficient. J Comput Phys 7:120–133
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, Patton EG (2012) Momentum and scalar transport within a vegetation canopy following
atmospheric stability and seasonal canopy changes: the CHATS experiment. Atmos Chem Phys
12:5913–5935
Dwyer MJ, Patton EG, Shaw RH (1997) Turbulent kinetic energy budgets from a large-eddy simulationof airflow above and within a forest canopy. Boundary-Layer Meteorol 84:23–43
Finnigan J (2000) Turbulence in plant canopies. Annu Rev Fluid Mech 32:519–571
Finnigan JJ, Shaw RH, Patton EG (2009) Turbulence structure above a vegetation canopy. J Fluid Mech
637:387–424
Fitzmaurice L, Shaw RH, Paw U KT, Patton EG (2004) Three-dimensional scalar microfront systems ina large-eddy simulation of vegetation canopy flow. Boundary-Layer Meteorol 112:107–127
Gao W, Shaw RH, Paw U KT (1989) Observation of organized structure in turbulent flow within andabove a forest canopy. Boundary-Layer Meteorol 47:349–377
Gavrilov K, Morvan D, Accary G, Lyubimov MS (2013) Numrical simulation of coherent turbulent
structures and of passive scalar dispersion in a canopy sub-layer. Comput Fluids 78:54–62
Geier M, Schönherr M, Pasquali A, Krafczyk M (2015) The cumulant lattice Boltzmann equation inthree 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

13

70

T. Watanabe et al.

He X, Chen S, Doolen GD (1998) A novel thermal model for the lattice Boltzmann method in incompressible limit. J Comput Phys 146:282–300
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.
Boundary-Layer Meteorol 164:161–181
Iwata H, Harazono Y, Ueyama M, Sakabe A, Nagano H, Kosugi Y, Takahashi K, Kim Y (2015) Methane
exchange in a poorly-drained black spruce forest over permafrost observed using the eddy covariance technique. Agric For Meteorol 214–215:157–168
Katul GG, Albertson JD (1998) An investigation of higher-order closure models for a forested canopy.
Boundary-Layer Meteorol 89:47–74
Katul GG, Chang WH (1999) Principal length scales in second-order closure models for canopy turbulence. J Appl Meteorol 38(11):1631–1643
Katul G, Hsieh CI, Kuhn G, Ellsworth D (1997a) Turbulent eddy motion at the forest-atmosphere interface. J Geophys Res 102:13409–13421
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
Katul G, Kuhn G, Schieldge J, Hsieh CI (1997b) The ejection-sweep character of scalar fluxes in theunstable surface layer. Boundary-Layer Meteorol 83:1–26
Katul G, Peltola O, Grönholm T, Launiainen S, Mammarella I, Vesala T (2018) Ejective and sweeping
motions above a peatland and their role in relaxed-eddy-accumulation measurements and turbulent
transport modelling. Boundary-Layer Meteorol 169:163–184
Kobayashi H (2005) The subgrid-scale models based on coherent structures for rotating homogeneous turbulence and turbulent channel flow. Phys Fluids 17:045104
Kondo J, Watanabe T (1992) Studies on the bulk transfer coefficients over a vegetated surface with a multilayer energy budget model. J Atmos Sci 49:2183–2199
Krüger T, Kusumaatmaja H, Kuzmin A, Shardt O, Silva G, Viggen EM (2017) The lattice Boltzmann
method. Springer, Switzerland
Kuwagata T, Hamasaki T, Watanabe T (2008) Modeling water temperature in a rice paddy for agro-environmental research. Agric For Meteorol 148:1754–1766
Lallemand P, Luo LS (2003) Lattice Boltzmann method for moving boundaries. J Comput Phys
184:406–421
Maitani T, Seo T (1985) Turbulent transport of scalar quantities within and above a paddy field. BoundaryLayer Meteorol 33:197–208
Maitani T, Shaw RH (1990) Joint probability analysis of momentum and heat fluxes at a deciduous forest.
Boundary-Layer Meteorol 52:283–300
Mellor GL, Yamada T (1974) A hierarchy of turbulence closure models for planetary boundary layers. J
Atmos Sci 31:1791–1806
Meyers TP, Paw U KT (1986) Testing of a higher-order closure model for modelling airflow within andabove plant canopies. Boundary-Layer Meteorol 37:297–311
Meyers TP, Paw U KT (1987) Modelling the plant canopy micrometeorology with higher-order closure
principles. Agric For Meteorol 41:143–163
Murayama S, Takamura C, Yamamoto S, Saigusa N, Morimoto S, Kondo H, Nakazawa T, Aoki S, Usami T,
Kondo M (2010) Seasonal variations of atmospheric ­CO2, δ13C, and δ18O at a cool temperate deciduous forest in Japan: influence of Asian monsoon. J Geophys Res 115:D17304
Müller A, Miyazaki Y, Tachibana E, Kawamura K, Hiura T (2017) Evidence of a reduction in cloud condensation nuclei activity of water-soluble aerosols caused by biogenic emissions in a cool-temperate
forest. Sci Rep 7:8452
Nakagawa H, Nezu I (1977) Prediction of the contributions to the Reynolds stress from bursting events inopen-channel flows. J Fluid Mech 80:99–128
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
Onodera N, Idomura Y, Hasegawa Y, Nakayama H, Shimokawabe T, Aoki T (2021) Real-time tracer dispersion simulations in Oklahoma City using the locally mesh-refined lattice Boltzmann method. Boundary-Layer Meteorol 179:187–208
Poggi D, Katul GG, Vidakovic B (2011) The role of wake production on the scaling laws of scalar concentration fluctuation spectra inside dense canopies. Boundary-Layer Meteorol 139:83–95
Raupach MR (1981) Conditional statistics of Reynolds stress in rough-wall and smooth-wall turbulent
boundary layers. J Fluid Mech 108:363–382

13

Coherent Eddies Transporting Passive Scalars Through the Plant…

71

Raupach MR (1989) Applying Lagrangian fluid mechanics to infer scalar source distributions from concentration profiles in plant canopies. Agric For Meteorol 47:85–108
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, Thom AS (1981) Turbulence in and above plant canopies. Annu Rev Fluid Mech 13:97–129
Shan X (1997) Simulation of Rayleigh-Bénard convection using a lattice Boltzmann method. Phys Rev E
55:2780–2788
Shaw RH (1977) Secondary wind speed maxima inside plant canopy. J Appl Meteorol 16:514–521
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, Paw U KT, Gao W (1989) Detection of temperature ramps and flow structures at a deciduous forest site. Agric For Meteorol 47:123–138
Shaw RH, Tavanger J, Ward DP (1983) Structure of the Reynolds stress in a canopy layer. J Climate Appl
Meteorol 22:1922–1931
Shaw RH, Zhang XJ (1992) Evidence of pressure-forced turbulent flow in a forest. Boundary-Layer Meteorol 58:273–288
Siqueira M, Lai CT, Katul G (2000) Estimating scalar sources, sinks, and fluxes in a forest canopy using
Lagrangian, Eulerian, and hybrid inverse models. J Geophys Res 105:29475–29488
Su HB, Shaw RH, Paw U KT (2000) Two-point correlation analysis of neutrally stratified flow within andabove a forest from large-eddy simulation. Boundary-Layer Meteorol 94:423–460
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
Takata K, Emori S, Watanabe T (2003) Development of the minimal advanced treatments of surface interaction and runoff. Global Planet Change 38:209–222
Tanaka K (2001) Multi-layer model of CO2 exchange in a plant community coupled with the water budgetof leaf surfaces. Ecol Model 147:85–104
Toda M, Watanabe T, Sumida A, Yokozawa M, Hara T (2009) Foliage profiles of individual trees determine competition, self-thinning, biomass and NPP of a Cryptomeria japonica forest stand: a simulation
study based on a stand-scale process-based forest model. Ecol Model 220:2272–2280
Ueyama M, Takanashi S, Takahashi Y (2014) Inferring methane fluxes at a larch forest using Lagrangian,Eulerian, and hybrid inverse models. J Geophys Res Biogeosci 119:2018–2031
Yoshino M, Inamuro T (2003) Lattice Boltzmann simulations for flow and heat/mass transfer problems in athree-dimensional porous structure. Int J Numer Meth Fluids 43:183–198
Watanabe T (1993) The bulk transfer coefficients over a vegetated surface based on K-theory and a 2ndorder closure model. J Meterol Soc Jpn 71:33–42
Watanabe T (1994) Bulk parameterization for a vegetated surface and its application to a simulation of nocturnal drainage flow. Boundary-Layer Meteorol 70:13–35
Watanabe T (2004) Large-eddy simulation of coherent turbulence structures associated with scalar ramps
over plant 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
Watanabe T, Mizutani K (1996) Model study on micrometeorological aspects of rainfall interception over anevergreen broad-leaved forest. Agric For Meteorol 80:195–214
Watanabe T, Shimoyama K, Kawashima M, Mizoguchi Y, Inagaki A (2020) Large-eddy simulation of neutrally-stratified turbulent flow within and above plant canopy using the central-moments-based lattice
Boltzmann method. Boundary-Layer Meteorol 176:35–60
Watanabe T, Yokozawa M, Emori S, Takata K, Sumida A, Hara T (2004) Developing a multilayered integrated numerical model of surface physics–growing plants interaction (MINoSGI). Global Change
Biol 10:963–982
Wilson NR, Shaw RH (1977) A higher order closure model for canopy flow. J Appl Meteorol 16:1197–1205
Zhuang Y, Amiro BD (1994) Pressure fluctuations during coherent motions and their effects on the budgetsof turbulent kinetic energy and momentum flux within a forest canopy. J Appl Meteorol 33:704–711
Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps andinstitutional affiliations.

13