Coherent Eddies Transporting Passive Scalars Through the Plant Canopy Revealed by Large-Eddy Simulations Using the Lattice Boltzmann Method

Abstract

A double-distribution-function lattice Boltzmann model for large-eddy simulations of a passive scalar field in a neutrally stratified turbulent flow is described. In simulations of the scalar turbulence within and above a homogeneous plant canopy, the model’s performance is found to be comparable with that of a conventional large-eddy simulation model based on the Navier–Stokes equations and a scalar advection–diffusion equation in terms of the mean turbulence statistics, budgets of the second moments, power spectra, and spatial two-point 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 the canopy 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 vice versa, 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 well as the transport of scalars emitted from the underlying soils to the atmosphere.

1 Introduction

A plant canopy produces a peculiar internal environment through physical and biological processes. The resulting microclimate exerts an essential influence on the energy balance, 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/from the atmosphere above serve to renew the microclimate inside the canopy and convey the physical 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 and Thom 1981; Finnigan 2000; Brunet 2020). The coherent eddies, which have vertical scales 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 of the 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 (or K-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 lower-level 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 energy balance 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 nature and mechanism of the non-local transport phenomena is still lacking, it is important to obtain 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 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), less attention 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 underneath crop 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), emissions of 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 processes of collision, forcing, and streaming. The macroscopic variables (fluid density, velocity, and pressure) derived from low-order moments of the distribution function are mathematically shown 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 a large-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 the advection 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 the macroscopic 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 representation of 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 in a 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 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 $\mathbf{x}$ = (${x}_{1}$, ${x}_{2}$, ${x}_{3}$) = ($x$, $y$, $z$) and $\mathbf{u}$ = (${u}_{1}$, ${u}_{2}$, ${u}_{3}$) = ($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, and the 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 (described later), 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 were reported 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, is hereafter 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 $\Delta x$, timestep $\Delta t$, reference fluid density ${\rho }_{0}$, and reference scalar density ${s}_{0}$. Therefore, in the normalized equations, these reference scales read $\Delta x=\Delta t={\rho }_{0}={s}_{0}=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); a summary is repeated here. The basic framework of the model is the central-moment-based lattice Boltzmann method in a three-dimensional 27-velocity (D3Q27) lattice system. Here and 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 ${f}_{ijk}\left(\mathbf{x},t\right)$ of discrete microscopic velocities $\left(ic,jc,kc\right)$, where the indices $i,j,k\in \left\{-\mathrm{1, 0},1\right\}$ represent the direction of the discretized velocities and $c=\Delta x/\Delta t$ (= 1) denotes the lattice velocity unit.

The distribution function and the macroscopic flow variables are updated by the following procedure. At every timestep, the central moments of the distribution function ${\kappa }_{pqr}$, are calculated by

$${\kappa }_{pqr}\left({\mathbf{x}}\text{, t}\right)=\sum_{i=-1}^{1}\sum_{j=-1}^{1}\sum_{k=-1}^{1}{\left(ic-u\right)}^{p}{\left(jc-v\right)}^{q}{\left(kc-w\right)}^{r}{f}_{ijk}\left({\mathbf{x}}\text{, t}\right), $$

(1)

where the indices $p,q,r\in \left\{\mathrm{0, 1},2\right\}$ 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, the post-collisional central moment ${\kappa }_{pqr}^{*}$ is obtained in terms of the pre-collisional central moments, the macroscopic fluid density and velocity, a net external force $\mathbf{F}=\left({F}_{\mathrm{x}},\boldsymbol{ }{F}_{\mathrm{y}},\boldsymbol{ }{F}_{\mathbf{z}}\right)$ working on the fluid, and the collisional relaxation coefficients ${\omega }_{1},\dots , {\omega }_{10}$ as

$${\kappa }_{pqr}^{*}\left(\mathbf{x},t\right)=\mathrm{func}\left({\varvec{\upkappa}},\rho ,\mathbf{u},\mathbf{F},{\omega }_{1},\dots , {\omega }_{10}\right). $$

(2)

The functional form of Eq. 2 is described in Watanabe et al. (2020). The inverse transform of Eq. 1 is then used to calculate the post-collision distribution function ${f}_{ijk}^{*}$, which is translated to the adjacent grid points in the next streaming step as

$${f}_{ijk}\left(x+ic\Delta t, y+jc\Delta t,z+kc\Delta t,t+\Delta t\right)={f}_{ijk}^{*}\left(x, y,z,t\right). $$

(3)

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

$$\rho =\sum_{i=-1}^{1}\sum_{j=-1}^{1}\sum_{k=-1}^{1}{f}_{ijk}, $$

(4)

$$\rho u=\sum_{i=-1}^{1}\sum_{j=-1}^{1}\sum_{k=-1}^{1}{\left(ic\right)f}_{ijk}+\frac{{F}_{\mathrm{x}}}{2}, $$

(5a)

$$\rho v=\sum_{i=-1}^{1}\sum_{j=-1}^{1}\sum_{k=-1}^{1}\left(jc\right){f}_{ijk}+\frac{{F}_{\mathrm{y}}}{2}, $$

(5b)

$$\rho w=\sum_{i=-1}^{1}\sum_{j=-1}^{1}\sum_{k=-1}^{1}\left(kc\right){f}_{ijk}+\frac{{F}_{\mathrm{z}}}{2}, $$

(5c)

$$p=\frac{\rho }{3}, $$

(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

$${\tau }_{\alpha \beta }=-\rho {\upsilon }_{\mathrm{SGS}}\left(\frac{\partial {u}_{\alpha }}{\partial {x}_{\beta }}+\frac{\partial {u}_{\beta }}{\partial {x}_{\alpha }}\right)+\frac{2}{3}{k}_{\mathrm{SGS}}{\delta }_{\alpha \beta } , $$

(7)

where ${\tau }_{\alpha \beta }$ represents the SGS stress, the indices $\alpha$ and $\beta$ indicate the spatial coordinates, ${\upsilon }_{\mathrm{SGS}}$ denotes the subgrid eddy viscosity, and ${\delta }_{\alpha \beta }$ is the Kronecker delta. The SGS kinetic energy denoted by ${k}_{\mathrm{SGS}}$ is estimated from

$$k_{{{\text{SGS}}}} = \rho c_{{\text{k}}} \left( {u_{\alpha } - \overline{\overline{u}} _{\alpha } } \right)^{2} ,$$

(8)

where ${c}_{\mathrm{k}}$ (= 1) is a model constant and $\overline{\overline{u}}_{\alpha }$ is a filtered velocity component obtained by

$$ \overline{\overline{u}} _{\alpha } = \frac{1}{{12}}\left( {u_{\alpha }^{{\text{B}}} + u_{\alpha }^{{\text{S}}} + u_{\alpha }^{{\text{W}}} + 6u_{\alpha } + u_{\alpha }^{{\text{E}}} + u_{\alpha }^{{\text{N}}} + u_{\alpha }^{{\text{T}}} } \right), $$

(9)

in which ${u}_{\alpha }^{\mathrm{B}}, {u}_{\alpha }^{\mathrm{S}},\dots$ 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}_{\alpha \beta }$ as

$${\upsilon }_{\mathrm{SGS}}={C}_{\mathrm{SGS}}{\left(\Delta x\right)}^{2}\sqrt{2{S}_{\alpha \beta }{S}_{\alpha \beta }}, $$

(10)

where ${C}_{\mathrm{SGS}}$ is a parameter determined according to the coherent-structure Smagorinsky model of Kobayashi (2005). The evaluated value of ${\upsilon }_{\mathrm{SGS}}$ is reflected in the lattice Boltzmann algorithm through the relaxation coefficient ${\omega }_{1}$ for a group of second-order central moments (Watanabe et al. 2020) as

$$\upsilon +{\upsilon }_{\mathrm{SGS}}=\frac{1}{3}\left(\frac{1}{{\omega }_{1}}-\frac{1}{2}\right). $$

(11)

The other relaxation coefficients of the central-moment-based collision model are set as ${\omega }_{2}={\omega }_{6}={\omega }_{7}={\omega }_{8}={\omega }_{10}={\omega }_{1}$ and ${\omega }_{3}={\omega }_{4}={\omega }_{5}={\omega }_{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

$$\frac{\partial \rho }{\partial t}+\frac{\partial \rho {u}_{\alpha }}{\partial {x}_{\alpha }}=0$$

(12)

and the filtered Navier–Stokes equations

$$\frac{\partial \rho {u}_{\alpha }}{\partial t}+\frac{\partial \rho {u}_{\alpha }{u}_{\beta }}{\partial {x}_{\beta }}=-\frac{\partial p}{\partial {x}_{\alpha }}+\frac{\partial }{\partial {x}_{\beta }}\left[\rho \left(\upsilon +{\upsilon }_{\mathrm{SGS}}\right)\left(\frac{\partial {u}_{\alpha }}{\partial {x}_{\beta }}+\frac{\partial {u}_{\beta }}{\partial {x}_{\alpha }}\right)-\frac{2}{3}{k}_{\mathrm{SGS}}{\delta }_{\alpha \beta }\right]+{F}_{\alpha }. $$

(13)

2.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. the Bhatnagar–Gross–Krook model. This is the most common and simplest collision model used in the lattice Boltzmann literature (Chen and Doolen 1998; Aidum and Clausen 2010). 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 scalar density in preliminary simulations using more efficient D3Q7 and D3Q15 lattices, while the most accurate D3Q27 lattice requires more memory resources.

In this framework, the spatio–temporal development of the distribution function of a scalar, ${g}_{l}$ ($l=0,\dots ,18$), is described by (Krüger et al. 2017)

$${g}_{l}\left(\mathbf{x}+{\mathbf{c}}_{l}\Delta t, t+\Delta t\right)={g}_{l}\left(\mathbf{x}, t\right)-\frac{1}{{\tau }_{g}}\left({g}_{l}\left(\mathbf{x}, t\right)+\frac{1}{2}{Q}_{l}\left(\mathbf{x}, t\right)-{g}_{l}^{eq}\left(\mathbf{x}, t\right)\right)+{Q}_{l}\left(\mathbf{x}, t\right), $$

(14)

where ${\mathbf{c}}_{l}$ is one of the 19 microscopic velocities, defined as ${\mathbf{c}}_{0}=$ (0, 0, 0), ${\mathbf{c}}_{\mathrm{1,2}}=$ ($\pm c$, 0, 0), ${\mathbf{c}}_{\mathrm{3,4}}=$ (0, $\pm c$, 0), ${\mathbf{c}}_{\mathrm{5,6}}=$ (0, 0, $\pm c$), ${\mathbf{c}}_{7-10}=$ ($\pm c$, $\pm c$, 0),${\mathbf{c}}_{11-14}=$ ($\pm c$, 0, $\pm c$), ${\mathbf{c}}_{15-18}=$ (0, $\pm c$, $\pm c$), ${Q}_{l}$ is the increment of ${g}_{l}$ attributable to the scalar source, and ${g}_{l}^{eq}$ 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 from the scalar source, the function ${g}_{l}$ gradually approaches ${g}_{l}^{eq}$ at a time-rate defined by the collisional relaxation time ${\tau }_{g}$ (Guo et al. 2002). The equilibrium distribution function is given by (Chopard et al. 2009)

$${g}_{l}^{eq}={E}_{l}\, s\left[1+3{\mathbf{c}}_{l}\bullet \mathbf{u}+\frac{9}{2}{\left({\mathbf{c}}_{l}\bullet \mathbf{u}\right)}^{2}-\frac{3}{2}\mathbf{u}\bullet \mathbf{u}\right], $$

(15)

where $s$ denotes the macroscopic scalar density and ${E}_{l}$ denotes the weighing constants ${E}_{0}=1/3$, ${E}_{1-6}=1/18$, and ${E}_{7-18}=1/36$. The source term in Eq. 14 is given by

$${Q}_{l}={E}_{l} \left[{q}_{\mathrm{s}}+3{q}_{\mathrm{s}}{\mathbf{c}}_{l}\bullet \mathbf{u}+3{r}_{s}{\mathbf{c}}_{l}\bullet {\mathbf{F}}_{\mathrm{t}}\right], $$

(16)

where ${q}_{\mathrm{s}}$ denotes the scalar source (or sink if negative), ${r}_{s}$ is the scalar mixing ratio defined as ${r}_{s}=s/\rho$, and ${\mathbf{F}}_{\mathrm{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=\sum_{l=0}^{18}{g}_{l}+\frac{1}{2}{q}_{\mathrm{s}}. $$

(17)

By considering one-half of the source term here and in the collision term in Eq. 14, the overall numerical accuracy is maintained at the second order both in time and space (Guo et al. 2002; Krüger et al. 2017).

The scalar diffusion resulting from unresolved eddies is incorporated via an SGS gradient–diffusion model as

$${\tau }_{\mathrm{s}\alpha }=-{\kappa }_{\mathrm{SGS}}\frac{\partial s}{\partial {x}_{\alpha }}, $$

(18)

where ${\tau }_{\mathrm{s}\alpha }$ represents the SGS scalar flux and ${\kappa }_{\mathrm{SGS}}$ denotes the subgrid eddy diffusivity, which is parametrized using the subgrid eddy viscosity and the subgrid turbulent Schmidt number $Sc$ as

$${\kappa }_{\mathrm{SGS}}= {\upsilon }_{\mathrm{SGS}}/Sc. $$

(19)

The effect of the SGS scalar flux is reflected in the simulations by including ${\kappa }_{\mathrm{SGS}}$ in the definition of the collisional relaxation time as

$$\kappa +{\kappa }_{\mathrm{SGS}}=\frac{1}{3}\left({\tau }_{g}-\frac{1}{2}\right), $$

(20)

where $\kappa$ 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

$$\frac{\partial s}{\partial t}+\frac{\partial {u}_{\alpha }s}{\partial {x}_{\alpha }}=\frac{\partial }{\partial {x}_{\alpha }}\left[\left(\kappa +{\kappa }_{\mathrm{SGS}}\right)\frac{\partial s}{\partial {x}_{\alpha }}\right]+{q}_{\mathrm{s}}. $$

(21)

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 computational domain is driven by a streamwise external force ${F}^{\mathrm{ext}}$, which is fixed at ${F}^{\mathrm{ext}}/\left({\rho }_{0}\Delta x/{\left(\Delta t\right)}^{2}\right)={10}^{-7}$. The top boundary is a rigid free-slip wall at which the scalar density is fixed to ${s}^{\mathrm{top}}/{s}_{0}=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 the flow and the rigid bottom boundaries of the flow and scalar fields are represented by means of halfway specular reflection, in which distribution functions that stream into the domain from the wall are given by

$${f}_{ij\overline{k} }\left(x,y,{z}_{\mathrm{b}},t+\Delta t\right)={f}_{ij{k}^{{{\prime}}}}^{*}\left(x-ic\Delta t,y-jc\Delta t,{z}_{\mathrm{b}},t\right), $$

(22a)

$${g}_{\overline{l} }\left(x,y,{z}_{\mathrm{b}},t+\Delta t\right)={g}_{{l}^{{{\prime}}}}^{*}\left(x-{c}_{{l}^{{{\prime}}}x}\Delta t,y-{c}_{{l}^{{{\prime}}}y}\Delta t,{z}_{\mathrm{b}},t\right), $$

(22b)

where ${z}_{\mathrm{b}}$ is the height of the highest or lowest fluid node, and $\overline{k }$ and $\overline{l }$ denote inwards directions opposite to the directions towards the wall (${k}^{^{\prime}}$ and ${l}^{^{\prime}}$). 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$, $\rho$, $p$, $s$) are collocated at each grid node. 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 plant elements on the flow and scalar fields are represented by an instantaneous drag force ${\mathbf{F}}_{\mathrm{d}}$ and a source of the scalar ${q}_{\mathrm{s}}$, parametrized as

$${\mathbf{F}}_{\mathrm{d}}=-\rho {c}_{\mathrm{d}}a\left|\mathbf{u}\right|\mathbf{u}, $$

(23a)

$${q}_{\mathrm{s}}=-{c}_{\mathrm{s}}a\left|\mathbf{u}\right|\left(s-{s}_{\mathrm{c}}\right), $$

(23b)

where ${c}_{\mathrm{d}}$ and ${c}_{\mathrm{s}}$ are, respectively, the drag coefficient and the scalar exchange coefficient of canopy elements, $a$ is the leaf area density, and ${s}_{\mathrm{c}}$ is the scalar density at leaf surfaces. The canopy parameters are given as $ah=2$ (i.e. the leaf area index LAI = 2), ${c}_{\mathrm{d}}=0.2$, ${c}_{\mathrm{s}}/{c}_{\mathrm{d}}=0.2$, and ${s}_{\mathrm{c}}=0$ to enable idealized simulations excluding any complexity arisen from the canopy morphology and plant physiological factors. The frictional force and the scalar source from the underlying ground surface are given by

$${\mathbf{F}}_{\mathrm{g}}=-{\rho }_{\mathrm{b}}{C}_{\mathrm{Mg}}\left|{\mathbf{u}}_{\mathrm{b}}\right|{\mathbf{u}}_{\mathrm{b}}, $$

(24a)

$${q}_{\mathrm{sg}}=-{C}_{\mathrm{sg}}\left|{\mathbf{u}}_{\mathrm{b}}\right|\left({s}_{\mathrm{b}}-{s}_{\mathrm{g}}\right), $$

(24b)

where ${C}_{\mathrm{Mg}}$ and ${C}_{\mathrm{sg}}$ are the bulk transfer coefficients of momentum and the scalar, respectively, ${s}_{\mathrm{g}}$ 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 (${z}_{0\mathrm{g}}$) and the scalar (${z}_{\mathrm{sg}}$), which are prescribed as ${z}_{0\mathrm{g}}/h={10}^{-3}$ and ${z}_{\mathrm{sg}}/{z}_{0\mathrm{g}}=0.1$, respectively, throughout the bottom boundary. The total force in Eq. 16 is thus given as ${\mathbf{F}}_{\mathrm{t}}={\mathbf{F}}^{\mathrm{ext}}+{\mathbf{F}}_{\mathrm{d}}+{\mathbf{F}}_{\mathrm{g}}-\nabla p$. The macroscopic velocity satisfying Eqs. 5 and 23a is evaluated as (Guo and Zhao 2002)

$$\mathbf{u}=\frac{\stackrel{\sim }{\mathbf{u}}}{\frac{1}{2}+\sqrt{\frac{1}{4}+\frac{1}{2}{c}_{\mathrm{d}}a|\stackrel{\sim }{\mathbf{u}}|}}, $$

(25)

where $\stackrel{\sim }{\mathbf{u}}$ denotes the macroscopic velocity calculated without the drag force from Eq. 5, while the macroscopic scalar density satisfying Eqs. 17 and 23b is calculated from

$$s=\frac{\tilde{s }}{1+\frac{1}{2}{c}_{\mathrm{s}}a|\mathbf{u}|}, $$

(26)

in which $\tilde{s }$ denotes the scalar density calculated from Eq. 17 without the scalar source term. At the lowest nodes, the expressions ${c}_{\mathrm{d}}a$ and ${c}_{\mathrm{s}}a$ in these equations are replaced by $\left({c}_{\mathrm{d}}a+{C}_{\mathrm{Mg}}\right)$ and $\left({c}_{\mathrm{s}}a+{C}_{\mathrm{sg}}\right)$, respectively, to account for the friction and the scalar source from 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 from the Navier–Stokes model in the region well above the canopy.

The simulation is initialized using a uniform streamwise velocity component with a 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 10³ timesteps ($\approx 0.26h/{u}_{*}$) during the time periods between ± 10⁴ steps relative to the reference times, which are set every 10⁵ steps in the first 3 × 10⁶ 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 ($\stackrel{-}{}$) denotes the mean value, while a prime ($^{{\prime}}$) indicates deviations from the mean.

3 Results and Discussion

The results are presented in dimensionless form using the friction velocity ${u}_{*}$, which is defined at the canopy top by (Watanabe et al. 2020)

$${u}_{*}^{2}={F}^{\mathrm{ext}}\left({z}_{\mathrm{top}}-h\right)/{\rho }_{0}, $$

(27)

and the characteristic scale of the scalar ${s}_{*}$, estimated from

$${s}_{*}=-{\left(\overline{{s }^{{{\prime}}}{w}^{{{\prime}}}}+\overline{{\tau }_{\mathrm{s}3}}\right)}_{\mathrm{a}}/{u}_{*}, $$

(28)

where the subscript a denotes a mean value averaged over the height range $2\le z/h\le 7$.

3.1 Comparison with the Navier–Stokes Model

The performance of the lattice Boltzmann method is examined by comparing the results with those from a reference simulation by the established Navier–Stokes model. As the validation for the airflow field has been described elsewhere (Watanabe et al. 2020), the focus 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 the vertical scalar flux, calculated from the lattice Boltzmann and Navier–Stokes results. As a constant scalar density is imposed at the top boundary, the mean density and the variance from 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 an appreciable 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 top of 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, and the 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 monotonically with 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 the scalar 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.

Fig. 1

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

3.1.2 Budgets of the Scalar Variance and Flux

The budget equation for the resolved scalar variance is written as

$$\frac{\partial }{\partial t}\left(\frac{\overline{{s}^{{{\prime}}2}}}{2}\right)=0=-\overline{{s }^{{{\prime}}}{w}^{{{\prime}}}}\frac{\partial \overline{s}}{\partial z }-\frac{\partial \overline{{{s }^{{{\prime}}2}w}^{{{\prime}}}}}{\partial z}-\frac{\partial \overline{{s }^{{{\prime}}}{\tau }_{\mathrm{s}3}}}{\partial z}-{c}_{\mathrm{s}}a \overline{\left|\mathbf{u} \right| s {s}^{^{{\prime}}}}-\overline{{\varepsilon }_{\mathrm{s}}}-\overline{{G}_{\mathrm{div}}}, $$

(29)

where ${\varepsilon }_{\mathrm{s}}$ denotes the rate at which the resolved variance is converted to the SGS variance, which is ultimately dissipated by the molecular diffusion, and ${G}_{\mathrm{div}}$ is a term that arises from the compressibility of the fluid, written as

$${G}_{\mathrm{div}}=\frac{1}{2}\frac{\partial \overline{w}{s}^{{{\prime}}2}}{\partial z}+\frac{s{s}^{{{\prime}}}}{2}\frac{\partial {u}_{\alpha }}{\partial {x}_{\alpha }}.$$

(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 ${\varepsilon }_{\mathrm{s}}=-{{\tau }^{^{\prime}}}_{\mathrm{s}\alpha }\left(\partial {s}^{^{\prime}}/\partial {x}_{\alpha }\right)$, whereas the sum of $\overline{{\varepsilon }_{\mathrm{s}}}+\overline{{G}_{\mathrm{div}}}$ 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

$$\frac{\partial \overline{{s}^{{{\prime}}}{w}^{{{\prime}}}}}{\partial t}=0=-\overline{{w }^{{{\prime}}2}}\frac{\partial \overline{s}}{\partial z }-\frac{\partial \overline{{s }^{{{\prime}}}{w}^{{{\prime}}2}}}{\partial z}-{D}_{\mathrm{c}}^{sw}-{D}_{\mathrm{p}}^{sw}, $$

(31)

where the first two terms on the r.h.s. represent the flux production and the vertical transport by resolved turbulence, respectively. The term ${D}_{\mathrm{c}}^{sw}$ denotes the decorrelation attributable to the canopy drag and sink

$${D}_{\mathrm{c}}^{sw}={c}_{d}a\overline{\left| \mathbf{u} \right|ws^{{\prime}}}+{c}_{\mathrm{s}}a \overline{\left|\mathbf{u} \right|s{w}^{{\prime}}}, $$

(32)

and the term ${D}_{\mathrm{p}}^{sw}$, which is evaluated as the residual of Eq. 31 for both models, represents the other decorrelation effects caused by pressure and the SGS components.

Figure 2 shows profiles of each term in these budget equations calculated from the results 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 the variance budget in redistributing the flux from above the canopy to the interior, resulting in a slower equilibration (attained at $z/h\gtrsim 4$) than the scalar variance. The production terms in both budgets peak sharply at the canopy top, but decrease quickly with depth 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 variance and 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 the parametrizations 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 constant sink/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 second-order moments of scalar turbulence in plant canopies.

Fig. 2

Profiles of the terms in the budget equations of a scalar variance and b scalar flux simulated by the lattice 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 at selected 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 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 the largest wavenumbers. Although the horizontal derivatives are calculated by the pseudo-spectral 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 in the regions of sharp scalar fronts caused by intense vertical motions. Adopting a less-dispersive 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$\Delta$; this considerably enhances the numerical stability, but also requires a finer grid system than 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.06 $h/{\lambda }_{y}$, whereas the streamwise spectra gradually increase with decreasing wavenumber, indicating streamwise structures in the periodic computational domain.

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; ${L}_{x}$ and ${L}_{y}$ are the streamwise and lateral dimensions of the computational domain; ${\lambda }_{x}$ and ${\lambda }_{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)

3.1.4 Two-point Autocorrelation Functions of the Scalar Density

The spatial structure of the coherent eddies that induce energetic perturbations at a given reference height ${z}_{r}$ 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

$$ R_{{ab}} \left( {r_{x} ,r_{y} ,z,z_{r} } \right) = \frac{{\overline{{a^{\prime } \left( {x + r_{x} ,{\text{~}}y + r_{y} ,z} \right){\text{~}}b^{\prime } \left( {x,{\text{~}}y,z_{r} } \right)}} }}{{\left[ {\overline{{a^{{\prime 2}}\left( z \right)}} \, \overline{{b^{{\prime 2}}\left( {z_{r} } \right)}} } \right]^{{1/2}} }} $$

(33)

where $a$ and $b$ denote the variables, and ${r}_{x}$ and ${r}_{y}$ denote the streamwise and lateral separation, respectively. Figure 4 shows the spatial distribution of the scalar’s autocorrelation function ${R}_{\mathrm{ss}}$ for a reference height just above the canopy at ${z}_{r}/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 narrow in 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, inclined core region residing near the canopy top, a peripheral positive correlation with ${R}_{\mathrm{ss}}>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 streamwise-elongated and downstream-tilted features of ${R}_{\mathrm{ss}}$ are similar to those of ${R}_{uu}$ 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.

Fig. 4

Spatial distribution of two-point autocorrelation functions of the scalar density, ${R}_{ss}$, 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

The above comparisons confirm that the lattice Boltzmann method is as capable as the Navier–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, the characteristics 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 ${R}_{u\mathrm{s}}$ exhibits a very large, streamwise-tilted pattern, and that of ${R}_{w\mathrm{s}}$ is vertically coherent and horizontally compact. These patterns of ${R}_{u\mathrm{s}}$ and ${R}_{w\mathrm{s}}$ closely resemble the spatial structures of ${R}_{uu}$ and ${R}_{ww}$ reported by Watanabe et al. (2020) and indicate that the 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}^{^{\prime}}>0$) caused by eddies with ${u}^{^{\prime}}>0$ and ${w}^{^{\prime}}<0$ (Fig. 5c). As demonstrated by Watanabe (2004), when a high-speed downward (sweep) motion carrying air of high scalar concentration from 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 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 is decelerated 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, a combination 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 long tail of the same structure may cause a scalar-sweep event near the canopy floor, as will be discussed in Sect. 3.3.

Fig. 5

Vertical–streamwise distribution of two-point cross-correlation functions of a streamwise velocity component (${R}_{us}$) and b vertical velocity component (${R}_{ws}$), referenced to the scalar density at a height just above the canopy ($z/h$ = 1.05). In c, these functions are also displayed as vectors overlaid on the map of the two-point autocorrelation function of scalar density shown in Fig. 4c. The interval of line contours is 0.1

3.2.2 Quadrant Analysis

The joint probability density function (p.d.f.) of ${s}^{^{\prime}}$ and ${w}^{^{\prime}}$ for four different heights near the canopy top is displayed in Fig. 6. The figure also shows bold-line contours of the third-order truncated Gram–Charlier distribution (Antonia and Atkinson 1973; Nakagawa and Nezu 1977), given by

Fig. 6

Joint p.d.f.’s of normalized fluctuations of the scalar density and the vertical velocity component at four 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

$$p\left(\widehat{s},\widehat{w}\right)=G\left(\widehat{s},\widehat{w}\right)\left[1+\sum_{p+q=3}{\left(-1\right)}^{p+q}\frac{{m}_{pq}}{p!q!}\frac{{\partial }^{p+q}G\left(\widehat{s},\widehat{w}\right)}{\partial {\widehat{s}}^{p}{\widehat{w}}^{q}}\right], $$

(34)

where $\widehat{s}={s}^{^{\prime}}/{\sigma }_{\mathrm{s}}$ and $\widehat{w}={w}^{^{\prime}}/{\sigma }_{w}$ are the variables normalized by the corresponding standard deviations and $({m}_{pq}= \overline{{\widehat{s}}^{p} {\widehat{w}}^{q}})$ ($p,q\in \left\{\mathrm{0,1},\mathrm{2,3}\right\}$; $p+q=3$) are the normalized third-order central moments, or equivalently the normalized third-order cumulants. The function $G\left(\widehat{s},\widehat{w}\right)$ represents the two-dimensional Gaussian

$$G\left(\widehat{s},\widehat{w}\right)=\frac{1}{2\pi \sqrt{1-{R}_{sw}^{2}}}\mathrm{exp}\left[-\frac{{\widehat{s}}^{2}-2{R}_{sw}\widehat{s}\widehat{w}+{\widehat{w}}^{2}}{2\left(1-{R}_{sw}^{2}\right)}\right], $$

(35)

with ${R}_{sw}=\overline{\widehat{s}\widehat{w}}$ denoting the correlation coefficient between ${s}^{^{\prime}}$ and ${w}^{^{\prime}}$. As the covariance $\overline{{s}^{^{\prime}}{w}^{^{\prime}}}$ is negative in the present simulation, each quadrant of the coordinate is classified as follows:

Quadrant I ($\widehat{s}>0$, $\widehat{w}>0$) Outward interaction,

Quadrant II ($\widehat{s}<0$, $\widehat{w}>0$) Ejection,

Quadrant III ($\widehat{s}<0$, $\widehat{w}<0$) Inward interaction,

Quadrant IV ($\widehat{s}>0$, $\widehat{w}<0$) Sweep.

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 well-known fact that sweeps dominate over ejections in the vertical turbulent transport at these heights 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 ($0<\widehat{w}<1$) spreads in the negative region of the abscissa towards the minimum of the scalar perturbation at each level. A similar tendency has also been observed for the joint p.d.f. of temperature and ${w}^{^{\prime}}$ measured in a deciduous forest (Maitani and Shaw 1990), as well as that of ${u}^{^{\prime}}$ and ${w}^{^{\prime}}$ 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 to the 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}^{{\prime}}$ and ${w}^{^{\prime}}$ in relation to the eddy motions in each quadrant (Raupach 1981; Katul et al. 1997a, b, 2018; Cava et al. 2006). From this idealized 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 in the second moments calculated for each quadrant are analytically related to the third-order moments by (Nakagawa and Nezu 1977; Raupach 1981)

$${\overline{\left(\widehat{s}\widehat{w}\right)}}_{\mathrm{I}}-{\overline{\left(\widehat{s}\widehat{w}\right)}}_{\mathrm{III}}=\frac{1}{6\sqrt{2\pi }}\left[3\left({m}_{21}+{m}_{12}\right)-{R}_{sw}\left({m}_{30}+{m}_{03}\right)\right], $$

(36a)

$${\overline{\left(\widehat{s}\widehat{w}\right)}}_{\mathrm{IV}}-{\overline{\left(\widehat{s}\widehat{w}\right)}}_{\mathrm{II}}=\frac{1}{6\sqrt{2\pi }}\left[3\left({m}_{21}-{m}_{12}\right)-{R}_{sw}\left({m}_{30}-{m}_{03}\right)\right], $$

(36b)

$${\overline{\left(\widehat{s}\widehat{s}\right)}}_{\mathrm{I}}-{\overline{\left(\widehat{s}\widehat{s}\right)}}_{\mathrm{II}}-{\overline{\left(\widehat{s}\widehat{s}\right)}}_{\mathrm{III}}+{\overline{\left(\widehat{s}\widehat{s}\right)}}_{\mathrm{IV}}=\frac{2}{3\sqrt{2\pi }}{m}_{30}, $$

(37a)

$${\overline{\left(\widehat{w}\widehat{w}\right)}}_{\mathrm{I}}+{\overline{\left(\widehat{w}\widehat{w}\right)}}_{\mathrm{II}}-{\overline{\left(\widehat{w}\widehat{w}\right)}}_{\mathrm{III}}-{\overline{\left(\widehat{w}\widehat{w}\right)}}_{\mathrm{IV}}=\frac{2}{3\sqrt{2\pi }}{m}_{03}. $$

(37b)

The third-order mixed moments, ${m}_{12}$ and ${m}_{21}$, which are related to the vertical transport of the variance $(\overline{{s}^{{\prime}2}})$ and the covariance $\overline{{s}^{^{\prime}}{w}^{^{\prime}}}$, respectively, are thus written in terms of the second moments in each quadrant as

$${m}_{21}=\sqrt{2\pi }\left[{\overline{\left(\widehat{s}\widehat{w}\right)}}_{\mathrm{I}}-{\overline{\left(\widehat{s}\widehat{w}\right)}}_{\mathrm{II}}-{\overline{\left(\widehat{s}\widehat{w}\right)}}_{\mathrm{III}}+{\overline{\left(\widehat{s}\widehat{w}\right)}}_{\mathrm{IV}}\right]+\frac{\sqrt{2\pi }}{2}{R}_{sw}\left[{\overline{\left(\widehat{s}\widehat{s}\right)}}_{\mathrm{I}}-{\overline{\left(\widehat{s}\widehat{s}\right)}}_{\mathrm{II}}-{\overline{\left(\widehat{s}\widehat{s}\right)}}_{\mathrm{III}}+{\overline{\left(\widehat{s}\widehat{s}\right)}}_{\mathrm{IV}}\right], $$

(38a)

$${m}_{12}=\sqrt{2\pi }\left[{\overline{\left(\widehat{s}\widehat{w}\right)}}_{\mathrm{I}}+{\overline{\left(\widehat{s}\widehat{w}\right)}}_{\mathrm{II}}-{\overline{\left(\widehat{s}\widehat{w}\right)}}_{\mathrm{III}}-{\overline{\left(\widehat{s}\widehat{w}\right)}}_{\mathrm{IV}}\right]+\frac{\sqrt{2\pi }}{2}{R}_{sw}\left[{\overline{\left(\widehat{w}\widehat{w}\right)}}_{\mathrm{I}}+{\overline{\left(\widehat{w}\widehat{w}\right)}}_{\mathrm{II}}-{\overline{\left(\widehat{w}\widehat{w}\right)}}_{\mathrm{III}}-{\overline{\left(\widehat{w}\widehat{w}\right)}}_{\mathrm{IV}}\right]. $$

(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. 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 effect of skewness), the moments ${m}_{12}$ (or ${m}_{21}$) essentially represent the difference between the contributions of eddy motions with positive and negative perturbations of ${w}^{^{\prime}}$ (or ${s}^{^{\prime}}$) to the covariance. Thus, non-Gaussian motions such as coherent eddies developing in the canopy turbulence, which cause asymmetry in the variances and covariances associated with downward and upward motions (especially sweeps and ejections), is responsible for the large contributions of the turbulent transport terms to the budgets of $\overline{{s}^{{\prime}2}}$ and $\overline{{s}^{^{\prime}}{w}^{^{\prime}}}$ (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 and scalars 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}$ and ${m}_{12}$ are directly related to the imbalance between sweeps and ejections ${\overline{\left(\widehat{s}\widehat{w}\right)}}_{\mathrm{IV}}-{\overline{\left(\widehat{s}\widehat{w}\right)}}_{\mathrm{II}}$ when the term ${R}_{sw}\left({m}_{30}-{m}_{03}\right)$ is neglected and a linear relationship between ${m}_{21}$ and ${m}_{12}$ 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 is because the term ${R}_{sw}\left({m}_{30}-{m}_{03}\right)$ is not negligible—${m}_{30}$ 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 ${m}_{12}/{m}_{21}$ to vary with height (Fig. 7b).

Fig. 7

Vertical profiles of a the normalized third moments of scalar density and vertical velocity components and b the ratio of ${m}_{12}$ to ${m}_{21}$. Marks denote the values calculated directly from the lattice Boltzmann results, and lines represent the estimations from Eqs. 37 and 38

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, which cause 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 the scalar 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 vertical motions 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 variance 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 numerical simulation 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.

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^{{\prime}}/{u}_{*}$, $v^{{\prime}}/{u}_{*}$), and colour shading indicates the vertical velocity component $w^{{\prime}}/{u}_{*}$ in (a) and the scalar fluxes $s^{{\prime}}w^{{\prime}}/\left({u}_{*}{s}_{*}\right)$ by sweep (blue) and ejection (red) motions in (b) with the sign of the ejection flux reversed

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

$$\langle c^{{\prime}}\rangle \left({r}_{x},{r}_{y},z,{z}_{r},\tau \right)=\frac{\sum_{j=1}^{n}\sum_{i=1}^{{m}_{j}}{c}^{{{\prime}}}\left({x}_{r}^{i}+{r}_{x},{y}_{r}^{i}+{r}_{y}, z,{t}_{j}+\tau \right)}{\sum_{j=1}^{n}{m}_{j}}, $$

(39)

where $\langle c^{{\prime}}\rangle$ denotes an ensemble average of any variable $c^{{\prime}}$ sampled at time ${t}_{j}+\tau$ (with $\tau$ representing a time lag) around detection positions $({x}_{r}^{i},{y}_{r}^{i},{z}_{r})$ ($i=1,\dots ,{m}_{j}$), where the detection conditions are satisfied in the j-th sample obtained at time ${t}_{j}$, ${r}_{x}$ and ${r}_{y}$ are the streamwise and lateral displacements from the detection position, respectively, $n$ is the number of samples, and ${m}_{j}$ is the number of detected events in the j-th sample. The height of the detection position is set close to the canopy floor at ${z}_{r}/h=0.15$, and the conditions for the sweep and ejection events are set as follows:

Scalar sweep: $w^{{\prime}}<0$ and ${s}^{^{\prime}}w^{{\prime}}/\left({\sigma }_{\mathrm{s}}{\sigma }_{w}\right)<-2$,

Scalar ejection: $w^{{\prime}}>0$ and ${s}^{^{\prime}}w^{{\prime}}/\left({\sigma }_{\mathrm{s}}{\sigma }_{w}\right)<-2$.

Figures 9 and 10 show vertical slices of the no-lag ($\tau =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 of the 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 autocorrelation of ${s}^{^{\prime}}$ 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.

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 (${z}_{r}/h=0.15$): a streamwise velocity component $\langle u^{{\prime}}\rangle /{u}_{*}$, b vertical velocity component $\langle w^{{\prime}}\rangle /{u}_{*}$, c pressure $\langle p^{{\prime}}\rangle /\left(\rho {u}_{*}^{2}\right)$, and d scalar density $\langle {s}^{^{\prime}}\rangle /{s}_{*}$ with vectors of ($\langle u^{{\prime}}\rangle /{u}_{*}$, $\langle w^{{\prime}}\rangle /{u}_{*}$). Contour intervals are a 0.2, b 0.1, c 0.2, and d 0.3

Fig. 10

Same as Fig. 9, but for scalar ejection events

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 is a 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 ${r}_{x}/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 vertically coherent, 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 the detection position, while that of the downdraft appears downstream and is slightly lifted from 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 streamwise velocity component, lies above an area of relatively lower scalar density that is accompanied 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 the detection position, whereas the updraft centre in the ejection case is located exactly over the detection position (Fig. 11). The most active parts of the eddy structures are streamwise-elongated (for sweeps) or horizontally compact (for ejections), rather than having transversally 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.

Fig. 11

Horizontal cross-sections at a height just above the canopy ($\mathrm{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 (${z}_{r}/h=0.15$). Vectors represent the horizontal velocity components ($\langle {u}^{{{\prime}}}\rangle /{u}_{*}$, $\langle {v}^{{{\prime}}}\rangle /{u}_{*}$), and coloured shading indicates the vertical velocity component $\langle {w}^{{{\prime}}}\rangle /{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 of the 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.

Fig. 12

Vertical–streamwise cross-sections of time-lagged composite structures for eddies causing scalar sweep events at a height close to the canopy floor (${z}_{r}/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 $\langle {s}^{{{\prime}}}\rangle /{s}_{*}$ and vectors represent ($\langle u^{{\prime}}\rangle /{u}_{*}$, $\langle w^{\prime}\rangle /{u}_{*}$)

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 as 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 the canopy 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. $\langle {u}^{^{\prime}}\rangle \langle {w}^{{\prime}}\rangle <0$; Fig. 13c), the downdraft of the sweep 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 of ejected lower scalar density is then lifted further by the remaining updraft and stretched by the mean shear (Fig. 13e).

Fig. 13

Same as Fig. 12, but for scalar ejection events

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 and 10d). The maximum values of the downward velocity component and pressure are also 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 case are observed at the level of the canopy height (Fig. 10a, b). The vertical coherency of these variables is more pronounced than for the sweep case, in which the pattern of the vertical velocity 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 around the 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 ensemble-averaging 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. This similarity 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 originated from 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 the two-dimensional roller-like structures themselves remain to contribute the Reynolds stress near 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 the canopy floor or if additional mechanisms are involved.

4 Conclusions

We have described a double-distribution-function lattice Boltzmann model for simulating a passive 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 performance of 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 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 from the 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 the sweep 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 from within and above plant canopies.

The scalar sweep events near the canopy floor are a direct consequence of the penetration of 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 in-phase at levels from near the canopy floor to approximately twice the canopy height at their maximum strength, when scalar ejection events are induced near the canopy floor.

These coherent eddies contribute to the non-local transport of the scalars emitted/absorbed by 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 the joint 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 is unchanged, the updrafts generated near or above the canopy top play a more significant role in transporting those scalars that have their primary sources at the soil surface (e.g., respiratory CO₂, CH₄, 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 of taking 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 non-local transport of scalars.

The characteristics of the coherent eddies described here were deduced from a single case of 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 to thermally non-neutral cases with marginal modifications to the code and would be feasible for further investigations on the non-local transport phenomena in wider conditions.