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.