Realistic Forests and the Modeling of Forest‐Atmosphere Exchange

Forests cover nearly a third of the Earth's land area and exchange mass, momentum, and energy with the atmosphere. Most studies of these exchanges, particularly using numerical models, consider forests whose structure has been heavily simplified. In many landscapes, these simplifications are unrealistic. Inhomogeneous landscapes and unsteady weather conditions generate fluid dynamical features that cause observations to be inaccurately interpreted, biased, or over‐generalized. In Part I, we discuss experimental, theoretical, and numerical progress in the understanding of turbulent exchange over realistic forests. Scalar transport does not necessarily follow the flow in realistic settings, meaning scalar quantities are rarely at equilibrium around patchy forests, and significant scalar fluxes may form in the lee of forested hills. Gaps and patchiness generate significant spatial fluxes that current models and observations neglect. Atmospheric instability increases the distance over which fluxes adjust at forest edges. In deciduous forests, the effects of patchiness differ between seasons; counter intuitively, eddies reach further into leafy canopies (because they are rougher aerodynamically). Air parcel residence times are likely much lower in patchy forests than homogeneous ones, especially around edges. In Part II, we set out practical ways to make numerical models of forest‐atmosphere more realistic, including by accounting for reconfiguration and realistic canopy structure and beginning to include more chemical and physical processes in turbulence resolving models. Future challenges include: (a) customizing numerical models to real study sites, (b) connecting space and time scales, and (c) incorporating a greater range of weather conditions in numerical models.

Previous reviews of ecosystem-scale turbulent transport around forests have discussed land-atmosphere interactions more generally (e.g., Bou-Zeid et al., 2020;Fisher & Koven, 2020) and ecosystem exchange at long-term monitoring sites (e.g., Baldocchi, 2008;Hicks & Baldocchi, 2020;Oliphant, 2012). These reviews, and the work they discuss, attempt to estimate the net ecosystem exchange of carbon and other scalar quantities, sometimes for scaling up into global-system and other large-scale models. Other reviews focus on the detailed structure and statistics of the turbulence in forests and other vegetated landscapes. The early work in this area was motivated by applied problems in agriculture and forestry, such as minimizing crop lodging and the windthrow of trees (e.g., Inoue, 1955;Ruth & Yoder, 1953). In the 1970s and 1980s, pioneering investigations of the transport equations moved the field from empirical observations to obtaining precise knowledge about certain features of the turbulence, especially in horizontally uniform canopies on flat ground (Finnigan, 2000;Kaimal & Finnigan, 1994;Raupach & Thom, 1981;Raupach et al., 1996;Van Gardingen & Grace, 1991). Katul et al. (2013) supplement these reviews with their survey of scalar turbulence in uniform canopies. Lemone et al. (2019) and Brunet (2020) provide comprehensive background to boundary-layer measurements, including in plant canopies, and update earlier reviews of the flow in homogeneous plant canopies with results obtained using high-performance computers. Other recent reviews discuss the flow and transport of scalar quantities around forest edges (Belcher et al., 2012), the measurement and modeling of ecological processes such as evapotranspiration (e.g., Katul et al., 2012), and hill-induced dynamical effects on the flow and scalar fields around forests and other vegetation (e.g., Belcher et al., 2012;Finnigan et al., 2020).

The Scope of This Review
The majority of existing work considers forests whose structure has been "idealized" in some way. Often this idealization step occurs when physical and numerical models are configured, although occasionally study sites are chosen on the basis of strict criteria. While there is no single definition of "idealized" in the context of forest-atmosphere exchange, the majority of the literature assumes forests to be horizontally homogeneous, with closed canopies, no large patches or gaps, and without site specific morphological or ecological detail. Below, we describe a forest as "idealized" if most of these properties are assumed. Forests in the real world rarely fit these assumptions. For example, forests around the world are becoming increasingly fragmented (Bogaert et al., 2011;Fahrig, 2003;Taubert et al., 2018). Edge regions differ from the forest interior both in their mean local climate (e.g., less humid) and in the range of meteorological extremes they experience (Magnago et al., 2015). Even in intact ecosystems, most forest canopies comprise a patchwork of openings of many shapes and sizes, formed by senescence, disease, and windthrow (Hirons & Thomas, 2018; Whitmore, 1989). These gaps are significant ecologically and structurally. The floras of the northern temperate forests, for example, include many species that depend on gaps and patchiness (Fox, 1977;Tinya et al., 2009). The divergence between theory and reality causes, for example, estimates of ecosystem exchange to be inaccurate or biased (Aubinet et al., 2010;Baldocchi, 2008;Stoy et al., 2013;K. Wilson et al., 2002). This is no secret, and most studies end with calls for further experiments in a wider range of ecosystems, using more detailed techniques and models.
We focus our attention on ecosystem-scale exchange of momentum, energy, and mass between the atmosphere and heterogeneous forests. We split our review into two parts. Part I (Sections 2-4) discusses experimental, theoretical, and numerical progress in our understanding of turbulent exchange over realistic, patchy forests. Section 2 summarizes the methodology behind numerical investigations of forest-atmosphere exchange, and the turbulence structure around idealized forests, including the effects of atmospheric stability. The remainder of Part I discusses how structural heterogeneity and patchiness affects the exchange of momentum (Section 3) and scalar quantities (Section 4). We concentrate principally on forests on level ground, although Section 4 briefly discusses recent progress in investigating scalar transport around forested hills. (See Finnigan et al. [2020], for example, for a comprehensive treatment of terrain induced effects on land-atmosphere exchange.) Part II (Sections 5-6) sets out practical ways to make numerical models of forest-atmosphere more realistic, including by accounting for reconfiguration and realistic canopy structure (Section 5) and beginning to include more chemical and physical processes in turbulence resolving models (Section 6). We conclude in Section 7 and provide recommendations for further research under the broader goals of (a) customizing numerical models to real study sites, (b) connecting space and time scales, and (c) incorporating a greater range of weather conditions in numerical models.

Part I: Forest-atmosphere exchange around heterogeneous forests
We begin Part I by summarizing the important fluid dynamical properties of forest-atmosphere exchange in idealized landscapes before moving on to look at exchange in more realistic landscapes.
Forests are located in the atmospheric boundary layer (ABL), the layer of atmosphere that is directly influenced by the Earth's surface. Figure 1 presents a schematic of the structure of the daytime ABL. The top of the daytime ABL, at height z i , caps the mixed layer, within which variables such as humidity and potential temperature are approximately constant with height. The lowest 10% or so of the ABL is known as the atmospheric surface layer (ASL), analogous to the "inner region" in wall boundary layers. Above rough surfaces, such as forests and urban areas, the ASL can be further divided into the inertial sublayer (ISL) and the roughness sublayer (RSL) Figure 1. Sublayers of the daytime atmospheric boundary layer (ABL) over a forest. The figure follows the classification in Oke (1988) for an urban boundary layer. The height of the mixed layer, which typically accounts for around 90% of the daytime ABL height, is suppressed to aid presentation. The variables z i and h c denote the depth of the ABL and the mean height of the forest, respectively.
= 0, where P is the kinematic pressure, g is the gravitational acceleration, θ 0 is a reference temperature, θ v is the virtual potential temperature, and δ ij is the Kronecker delta (the second term on the right-hand side of Equation 2b is non-zero when i = 3). Capital letters denote the double-averaged quantities, which we refer to as the mean quantities-for example, the mean streamwise velocity component U. The term − ⟨̄̄⟩∕ is the dispersive flux of mean momentum, which accounts for spatial correlations in the time-averaged velocity field. The dispersive flux is usually assumed to be low in homogeneous forests and is therefore typically disregarded in numerical models (Patton & Finnigan, 2012). However, recent evidence suggests that the dispersive fluxes of momentum and scalar quantities can be significant around patchy forests (see Section 6.4). See Finnigan (2000) and Finnigan and Shaw (2008) for more formal discussion of the double-average method.

The Double-Average Method in Numerical Models
In applications as complex as forest-atmosphere exchange, numerical models can only ever serve as incomplete models, rather than as mini versions of the real world, or even as computational wind tunnels. Nonetheless, numerical models have emerged as one of the most powerful tools to investigate forest-atmosphere exchange, supplementing experimental observations, theory, and physical models. Their main advantage is that the flow is defined at every cell within the simulated domain, allowing time series of interest to be extracted easily. With post-processing and visualization software, numerical models provide a picture of the flow that is impossible to replicate with observations or physical models. Researchers have adopted various approaches to approximate solutions to the double-averaged equations, including first-order analytical closure (Finnigan et al., 2015, and references therein), modified gradient-diffusion theory (Zeng & Takahashi, 2000), Reynolds-averaged Navier-Stokes (RANS) solvers (Boudreault et al., 2015;Brunet, 2020, and references therein;Katul & Albertson, 1998), and large-eddy simulation (LES; Shaw & Schumann, 1992). RANS methods span first-order closure, "1.5-order" closure based on the transport of the turbulence kinetic energy and its (specific) rate of dissipation, and second-order and higher closure schemes . For a variety of applications dealing with land-atmosphere exchange, second-order closure appears to be sufficient . Direct numerical simulation (DNS) has recently been used for small, idealized plant canopies (Sharma & García-Mayoral, 2020a, 2020b. However, the computational expense of DNS means it can still be employed only for relatively low Reynolds number (Re) flow and that it remains unsuitable for ecosystem-scale investigations around forests. LES has emerged as the most popular method for investigating ecosystem-scale exchange around forests, although analytical closure schemes and RANS remain popular for situations that do not require the turbulence to be resolved.
Using LES, a low-pass filter of kernel G is applied to the flow field variables, where the overtilde denotes "resolved" field variables, whose prognostic equations are solved. The filter in Equation 3 directly resolves motion larger than the filter width Δ f , which is greater than or equal to the grid resolution Δ g , with ratios Δ f /Δ g = 1 or 2 common in practice (Basu & Porté-Agel, 2006;Geurts, 2003). The filtered equivalent of the averaged momentum equations in 2b are where ̃ accounts for forcings such as an imposed pressure gradient, the geostrophic wind, or distributed drag from vegetation and other obstacles smaller than the filter width, Δ f . The sub-grid scale (SGS) stress term τ ij accounts for unresolved scales, those smaller than Δ f , The SGS stresses and fluxes must be parametrized, for which a variety of schemes have been developed (Section 6.5; Gadde et al., 2021;Kirkil et al., 2012;Moser et al., 2021;Rozema et al., 2015). Once the LES model has been "spun up"-that is, the flow reaches a statistical equilibrium-the statistics of the flow variables are generated by averaging the resolved quantities, usually in space and time (e.g., Cai et al., 2008;B. Chen et al., 2019).
The turbulent perturbations are calculated as the departures from those averages (see Appendix A for an example). The main advantage of LES over RANS and analytical models is that LES expressly resolves the largest scales of turbulent motion in the simulated domain. Provided that Δ f and the model's forcings are chosen carefully, LES also allows a visualization and term-by-term analysis of the turbulent flow of air that cannot be achieved using observations or physical models. To avoid distracting switches between notation, we use the double-average notation throughout this manuscript. References to LES results and models should be interpreted as referring to the resolved variables, unless otherwise stated.

Representing Forests Through Distributed Drag
The aerodynamic drag of the forest (per unit mass of air) is accounted for through the term f i (m s −2 ) on the righthand side of Equation 2b. This term is the net sum of (a) the form drag from pressure differences either side of each plant element, and (b) the viscous boundary layers that develop over each plant element, where ν is the kinematic viscosity of air. This term is typically parametrized by spatially averaging the localized drag from the individual plant elements as where | | = ( ) 1∕2 and C d is a dimensionless drag coefficient. This parametrization is based on the aerodynamic drag equation, used in fluid mechanics and engineering applications. In the fluid mechanics literature, the drag equation (per unit mass) is usually written with a factor of 1/2, which originates from the formula for the kinetic energy of the fluid in front of the body. By meteorological convention, the factor of 1/2 is included in the drag coefficient and does not appear expressly here. This parametrization, which we refer to as the "distributed-drag parametrization," assumes the aerodynamic drag from the forest increases with the square of the velocity, as is the case around bluff bodies (Shaw & Schumann, 1992; N. R. Wilson & Shaw, 1977). The viscous component of the drag is usually neglected in the approximation of f i because form drag dominates in high-Re flow through forests (Shaw & Patton, 2003;Thom, 1971). The local forest density a(z) is typically assumed to be a function of the plant area density (PAD), the total one-sided plant area per unit layer volume (m 2 /m 3 ). The plant area index (PAI) is the PAD integrated over the height of the forest h c , that is, It is difficult to determine a robust approximation of the drag per unit volume, f i , in real forests. The parametrization is based on the form drag imposed by the frontal area of the various tree and plant elements (N. R. Wilson & Shaw, 1977). In closed canopies such as forests, it is typically assumed that more of the plant material is sheltered from the flow than in clumps of leaves or isolated trees, resulting a lower value of the product C d a(z) for closed canopies Sogachev & Panferov, 2006). Relatively low values of C d are often used in modeling studies to accommodate this shelter effect, with values ranging from 0.1 to 0.5 (Table A1). These values are usually assumed to be constant and invariant with height. Although sheltering effects, and therefore low values of C d , have been observed at moderate to high wind speeds (Amiro, 1990b;Koizumi et al., 2010), there is much more scatter at low wind speeds (Koizumi et al., 2010;Queck et al., 2012). This is probably because at low wind speeds, the leaves of the trees and plants are mostly relaxed, but will streamline temporarily in response to gusts, significantly altering the drag. At higher wind speeds, the plant elements spend more time in the streamlined position, so there is less scatter in the drag measurements during the passage of gusts. We discuss this reconfiguration in more detail in Section 5.

The Turbulence Structure Around Forests
As the air moves through a forest, momentum is transferred from the flow to the aerial parts of the plants and trees. This transfer reduces the streamwise wind speed throughout the depth of the canopy and a region of high wind shear forms around the crown top. The shear region is evidenced by an inflection in the mean streamwise windspeed profile, which is approximately exponential within the canopy and logarithmic above it (Finnigan, 2000;Raupach et al., 1996). Below the main crown, a secondary wind-speed maximum may occur (Shaw, 1977), especially near edges and in forests with sparse understories (Dupont et al., 2011). The source of this secondary maximum is not certain, although it may result from the turbulent transport of momentum from the upper canopy (Shaw, 1977) or mesoscale pressure gradients (Holland, 1989). The high shear around the crown top generates Kelvin-Helmholtz-type instabilities, which in turn generate coherent large eddies around the tops of the trees, analogous to the dominant processes in a plane mixing layer (Raupach et al., 1996). Using this mixing-layer analogy, Raupach et al. (1996) reduce canopy turbulence to a single length scale, = ℎ ∕( ∕ )ℎ , where ℎ is the mean streamwise velocity component U at z = h c . The shear length scale L s , which equates to around 0.5h c for medium density vegetation, provides a rough estimate of the diameter of the dominant turbulent eddies.
It is difficult to determine the value of L s exactly for real forests because dU/dz varies quickly with height around the crown top (i.e., d 2 U/dz 2 is not easily defined). However, L s can be estimated from the spectral energy peaks in wavelet cospectra of momentum and scalar quantities .
The second-order moments ⟨ ′2 ⟩ , ⟨ ′2 ⟩ , ⟨ ′ ′ ⟩, and the TKE increase with height within plant canopies, but are roughly constant above the canopy (Brunet, 2020;Raupach et al., 1996). The turbulent velocity components u′ and w′ are more correlated and skewed in the RSL than they are in the ISL above, where many velocity statistics display approximately Gaussian behavior. Streamwise and vertical skewness (Sk u and Sk w ) are approximately zero in the ISL but take values 0.5 ≤ Sk u ≤ 1 and −1 ≤ Sk w ≤ 0 in forests (Amiro, 1990a;Kruijt et al., 2000;Lee & Black, 1993;Raupach et al., 1996;Villani et al., 2003). The values of Sk u and Sk w may be higher, absolutely, at forest edges (Dupont & Brunet, 2008a). The association of positive Sk u and negative Sk w values indicates that turbulent transfer is dominated by strong but infrequent downward penetrations of air into the canopy, known as "sweep motions" (u′ > 0 and w′ < 0). Frequent upward motions of low-momentum air from within the canopy, known as "ejections" (u′ < 0 and w′ > 0), also account for a large proportion of the turbulent transfer (Shaw & Patton, 2003), although sweep motions contribute more to the total transfer of momentum around the crown top (Raupach et al., 1996;Shaw & Tavangar, 1983). Together sweep motions and ejections account for between 60% and 80% of the exchange of scalar quantities from homogeneous forest canopies to the atmosphere aloft (Gao et al., 1989). The greater magnitudes of the skewness statistics around forest edges reflect the significant differences in momentum and scalar exchange in patchy and gappy forests as compared to homogeneous canopies.
The mixing-layer analogy-where shear generated eddies around the crown top dominate turbulent exchangehas proved remarkably robust in forests and other vegetation. However, current understanding of canopy turbulence is far from complete. Further LES studies and targeted observations will help to reveal to what extent three-dimensional structures dominate the turbulence and the conditions under which the mixing-layer analogy breaks down when the canopy is patchy. Another knowledge gap concerns forest density. The vast majority of reported investigations have been on relatively dense tropical, Taiga, and northern temperate forests, usually in the growing season. In these forests, mixing-layer-type turbulence has been observed at low densities (Brunet, 2020). However, other sparse tree-dominated landscapes-say, the Guinea Savannah or temperate deciduous forests in winter-have similar mean densities (e.g., PAI 1-2; Xiao et al., 2017) but entirely different structural and meteorological drivers. These important landscapes remain poorly understood. For further discussion of plant-canopy turbulence see Finnigan (2000) for flow statistics and technical background, Finnigan et al. (2009) and Bailey and Stoll (2016) for the emergence of coherent fluid structures, Belcher et al. (2012) for scaling analysis in complex terrain, and Brunet (2020) for historical background and a review of recent studies in homogeneous plant canopies.

Atmospheric Stability and the Turbulence Structure
Air in the ASL is, on average, statically unstable during the day and stable at night, although the true dynamics are more complicated (Stull, 1991). The mixing-layer analogy of canopy turbulence assumes near-neutral conditions, with the dynamics controlled by the high shear in the mean wind velocity around the tops of the trees. However, in strongly stable or unstable conditions, the velocity shear can be much less influential (Brunet & Irvine, 2000;Lemone et al., 2019). As the ABL becomes more unstable, the turbulence structure around forests transitions from a shear-driven to a convection-driven regime, that is, thermal cells govern the flow dynamics, and the mixing-layer type turbulence becomes less prominent (Dupont & Patton, 2012;Lemone et al., 2019;Mahrt, 2000). Conversely, when the ABL is stable, the buoyancy of the air dampens vertical motion. Mixing-layer type coherent structures may still develop around a forest, but they are smaller and less energetic than those that form in near-neutral conditions (Dupont & Patton, 2012). ABL turbulence is generally more intermittent in space and time in stable conditions (Mahrt, 2014). In very stable conditions, fluxes of scalar quantities are driven by an interaction of mesoscale phenomena, such as gravity waves and nocturnal jets, and the local turbulence. This interaction is observed in the large, intermittent variations in temperature and CO 2 concentration around forest that occur when extended calm periods are interrupted by short bursts of intense turbulence (Aubinet, 2008;Heinesch et al., 2007;Wohlfahrt et al., 2005).
Simulations of flow around idealized forests in non-neutral conditions are beginning to add detail to this general picture. Patton et al. (2016) used LES to investigate the entire ABL over an interactive forest canopy, across five stability classes (near neutral to free convection). In strongly unstable conditions, thermal plumes may bubble up from the forest floor or the canopy top, and the vertical profiles of the atmosphere in the RSL approach those predicted by MOST. As stability decreases from near neutral to free convection, the dominant turbulent structures around the forest change from shear layer vortices to thermal plumes, and momentum and scalar fluxes become less correlated. It is not clear whether this transition is gradual, with shear and thermally driven structures coexisting, or sudden, with a flip between regimes upon reaching a stability threshold. Nascent research indicates the transition may be sudden (Brunet, 2020). In unstable conditions, scalar quantities are transported predominantly by thermal plumes. The Canopy Horizontal Array Turbulence Study (CHATS) observations provide the strongest observational evidence as to the effect of atmospheric stability on exchange in a tree dominated landscape (Dupont & Patton, 2012;Patton et al., 2011). These measurements found the velocity variances, momentum stresses, and momentum transport efficiency decrease as the atmosphere becomes less stable, indicating that thermal plumes become increasingly influential.
In stable conditions, wind shear at the crown top is higher than in unstable or neutral conditions, and momentum penetrates less deeply into the forest (Chaudhari et al., 2016;Nebenführ & Davidson, 2015;Su et al., 2008). Within the forest, turbulence is much weaker and more intermittent than in neutral or unstable conditions, and the pressure transport term of the TKE budget becomes more significant with increasing stability (Nebenführ & Davidson, 2015). In open forests, intermittent turbulence, driven by shear in the air aloft, may penetrate into the canopy and dramatically alter the distribution of scalars (Wharton et al., 2017). These intermittent events in stable conditions are not well understood and are seldom resolved well by numerical models because they do not result from resolved shear generated turbulence at the microscale or the MOST parametrizations used in larger scale models. The events are thought to result from nonturbulent "submeso" motions, which fall between the scale of the largest turbulent ABL eddies (∼O(100 m)) and the smallest γ-mesoscale (∼O(2 km); Mahrt, 2014). There is some evidence that parametrizations of submeso motions may be possible by analogy with "self-organized criticality," the tendency of dynamical systems to organize their microscopic behavior to be scale independent (Cava et al., 2019). However, this remains an active area of future research, requiring careful comparison between numerical model results and observations (Sun et al., 2015).
The atmospheric stability within forests often differs to that of the surrounding atmosphere. For example, on a sunny day, warm convective thermals may form above a forest, while the air within it remains stable (Ramos et al., 2004;Stull, 2006). At night, the situation is reversed when the forest crown loses heat through radiative cooling, forming a capping layer of very stable air around the tops of the trees that can decouple the forest air space from its surroundings (Nebenführ & Davidson, 2015;Paul-Limoges et al., 2017;X. Xu et al., 2015). The decoupling is sensitive to site specific meteorological conditions, such as local temperature gradients (Alekseychik et al., 2013;Russell et al., 2016). When the ABL and forest air space decouple, eddy-covariance measurements are difficult to interpret because CO 2 and other scalars can accumulate within the canopy, and the measured turbulent flux is unrepresentative of the forest as a whole (Aubinet, 2008;Massman & Lee, 2002). In hilly terrain, the decoupling can trigger drainage flows, transporting CO 2 , heat, and water vapor to and from forests (Finnigan et al., 2020;Sun et al., 2012;Wharton et al., 2017;X. Xu et al., 2015). Drainage flows and other below-canopy transport contribute a large proportion to the total scalar transport in certain meteorological conditions, such as during stable nights with weak winds (McHugh et al., 2017;Paul-Limoges et al., 2017).

Adaptation to Wind Loading
Wind affects all aerial parts of forests on time and space scales from those of the smallest eddy to those of climate variation and tree lifetimes. Many plants respond to wind by reducing their surface area to minimize drag, and therefore lower the likelihood of damage from wind loading. On longer timescales, plants may apportion biomass in such a way as to acclimate to the wind's prevailing direction and intensity, a process known as thigmomorphogenesis. Figure 2 presents examples of trees that show thigmomorphogenetic changes as a result of sustained wind loading. For example, trees at the upstream edge of a forest may have stiffer wood than those further inside the stand (Brüchert & Gardiner, 2006;Cucchi et al., 2004), and trees exposed to high winds develop lower crown densities, comprising a smaller number of smaller leaves (Telewski, 2009;Telewski & Pruyn, 1998). This adaptation is a trade-off for the trees, with increased resistance to wind loading and other environmental stresses coming at the expense of slower growth (Hirons & Thomas, 2018). For further background on the permanent response of trees to wind loading, see, for example, Telewski and Pruyn (1998), Telewski (2009), Hirons andThomas (2018).
Windbreaks, also known as shelter belts, exploit the ability of trees and shrubs to adapt to windy environments. Windbreaks have been used for centuries to shelter coastal buildings, to modify the microclimates around crops and grazing land, to protect woodland and wildlife, and, more recently, to modify pollutant transport around urban or industrial areas (Palmer et al., 1997;Pearce et al., 2021). Figure 3 shows a sketch of the streamlines around a windbreak for wind blowing left to right. The "quiet zone" is the region behind the barrier where the mean wind speeds are low and turbulence is dominated by recirculating eddies (Raine & Stevenson, 1977). The quiet zone extends x/h c ≈ 8-17 downstream of the windbreak, with higher values observed in conditions with low atmospheric turbulence. The "wake region" is above the quiet zone. The air in the wake region is highly turbulent because of the interaction of downward fluxes of momentum, and sharp velocity gradients as the flow adjusts around the barrier (Finnigan & Bradley, 1983). The sizes and energies of the vertical fluctuations are enhanced in the wake region compared to the incoming stream, whereas those of the vertical components are suppressed (McNaughton, 1989;H. Wang et al., 2001 and references therein).

Flow Adjustment Around Forest Edges
Compared to the flow around relatively isolated windbreaks, the flow around a forest is more complicated. Figure 4 presents a schematic of the statistical patterns in momentum transfer that emerge as the flow adjusts around a small forest in neutral atmospheric conditions (Belcher et al., 2003(Belcher et al., , 2012Dupont & Brunet, 2008a).
In the impact region, (i) in Figure 4, the forest acts as a step-change in porosity, inducing a pressure gradient to slow the flow. Just downstream of the forest edge is the adjustment region, (ii) in Figure 4, in which the drag from the trees decelerates the flow over a distance x A proportional to a canopy drag length scale L c = 1/C d a(z). The length scale L c emerges from quasi-inviscid solutions to the momentum equations in one-dimensional flow. It can be interpreted as a distance constant over which the velocity and drag adjust to balance the pressure gradient (Finnigan & Brunet, 1995). At the edge of homogeneous canopies, assuming constant shear and neutral atmospheric conditions, the canopy drag length scale L c is related to the shear length scale L s by where β = u * /U measured at the canopy top (z = h c ) (Brunet, 2020;Ross & Vosper, 2005). Because L c is inversely proportional to the plant density, the flow adjusts more quickly with increasing forest density, provided the density varies on scales greater than the volume averaging operation. (Planted forests are planted at a high density and then thinned repeatedly over the years to reach the final stocking density. It's not unusual for new forests to be planted at 1,000 stems per ha, or 3 m spacing (Forestry Commission, 2021; U.S. Department of Agriculture, 1982). Diameter at breast height measurements can then convert stand number density into basal area (i.e., m 2 (stem) ha −1 ).  Converting basal area into PAD requires knowledge of the allometry (i.e., size and shape) of the trees in the stand from standard models or remote sensing (Section 5.2, below). Variations from the planted-and-thinned density arise from introduction of rides for timber extraction, topographic features, and mortality from disease and windthrow.) The length scale L c is only an approximation for three-dimensional flow around real forests, for which the variables C d and a(z) may not be clearly defined (Section 5). Nonetheless, numerical simulations, field observations, and flume investigations of vegetation canopies show the adjustment distance x A downstream of a forest edge is indeed proportional to L c , with x A ≈ 4−6L c ≈ 8−10h c (Belcher et al., 2012;Morse et al., 2002;Rominger & Nepf, 2011;Yang, Raupach, et al., 2006). The canopy drag length scale permits an alternative expression for the PAI of a forest: PAI/h c = a(z) ≈ a, so that PAI ≈ h c /C d L c . This is only a rough approximation, because it assumes the plant area is evenly distributed in all directions. Modeling studies have used PAI values in the range 1-8, with most studies concentrating at the low-to-medium end of the range (Table A1). The relative scarcity of studies using PAI > 5 is surprising given that values in tropical and conifer forests often fall in the range 8-12 (Fleischbein et al., 2005;Lefsky et al., 1999). In the canopy-flow region, (iii) in Figure 4, the flow is fully adjusted to the presence of the forest. The canopy-shear layer (iv) is characterized by the shear-generated turbulent eddies that exchange most of the energy, mass and momentum between the forest and the atmosphere. These eddies are generated by processes analogous to those in a plane mixing layer (Section 2.5 above) (Raupach et al., 1996). Downstream of the adjustment region, an internal boundary layer may begin to develop above the trees, known as the roughness-change region (v).
If there is low vegetation or a large clearing in the lee of the forest, an exit region (vi) forms, within which mean wind speed increases, with a corresponding downwards flow, over a streamwise distance 1−4h c . The exit region may or may not contain a quiet zone of recirculating eddies, as behind windbreaks. Quiet zones are more likely to form in the lee of dense forests, especially those with an extensive understory, because the mean wind speed in the canopy flow region is typically low, meaning there is little exit flow at the leeward edge (Cassiani et al., 2008;Detto et al., 2008). However, because forest edges are typically more porous than purpose-built windbreaks, the recirculating eddies are more time intermittent than those behind windbreaks (Bergen, 1975;Cassiani et al., 2008). Recirculating eddies become increasingly common in forest lees as the forest density increases, although Cassiani et al. (2008) report a threshold PAI ≈ 6 after which further increases in density have little effect the turbulence structure in the lee. The formation and strength of the recirculating eddies appears to be relatively independent of a forest's foliage distribution (Ma et al., 2020), although the influence of edge morphology has not been carefully tested. A turbulent wake region (vii) occurs above the exit region. Because tree foliage creates substantial drag, the air above forests is typically more turbulent than that over the bare ground, sea, or fields typically found upstream of windbreaks. This higher turbulence decreases the distance over which the wake region reaches the ground behind forests as compared to windbreaks. The quiet zone is therefore relatively smaller behind forests (x/h c = 2-5) than behind windbreaks (x/h c = 8-17) and wake turbulence decays faster in the former. Banerjee and co-workers offer interesting arguments as to why the edge regions in Figure 4 appear to be so reproducible. Provided the forest is dense enough, the main features in Figure 4 are recovered from turbulently inviscid solutions of the momentum equations (Banerjee et al., 2013). This suggests that local forest structure and corresponding turbulence modifies the details of the features in Figure 4 but not their existence. In other words, the turbulently inviscid solution determines the spatial patterns in the mean flow field and turbulence displaces them locally or smooths them.

The Edge Regions in Context
Only about half of the world's remaining forest area, mostly in the Amazon and the Congo Basin, lies more than 500 m from the nearest edge (Haddad et al., 2015). In much of the Northern Hemisphere, forests are small and patchy because they are located close to areas where large populations of humans have lived for centuries. As an extreme example, approximately three-quarters of English woodland lies less than 100 m from the nearest edge (Riutta et al., 2014). Simple geometric considerations indicate that the area of forest subject to edge effects, the "edge region," is surprisingly large. Among two-dimensional shapes with the same area, a circle has the shortest perimeter, that is, it is the most compact shape in terms of its edge to area ratio. By approximating the plan of a forest stand as a circle, we obtain a lower limit to the area of the edge region from the annulus of width x A , where x A is the flow's adjustment distance. The lower limit for the ratio R 0 of the area of the edge region to the total stand area is therefore where A is the area of the forest stand and r is the radius of the equivalent-area circle. Taking R 0 > 1/2 in Equation 9 shows the edge region comprises over half the forest stand where Because most forest stands are not even approximately circular, the area subject to edge effects is substantially larger than this minimum. As a conservative heuristic for non-circular forests, we suggest edge effects dominate an area 25% greater than the area of the equivalent-area circle, therefore, where 1.25 × (6 + 4 √ 2) 2 ≈ 46 2 . Taking x A ≈ 8h c -the lower end of reported values of x A -provides a rule-ofthumb that edge effects dominate in forest stands where ∼ 3000ℎ 2 . Since forested areas are usually reported in hectares and canopy heights in meters, a dimensional version of the rule-of-thumb is For example, for a mature forest with a canopy height of 20 m, edge effects dominate for patches whose area are less than 120 ha. In many parts of the world (Haddad et al., 2015), edge effects dominate most of the forested area.

Patches and Gaps
Most real forests contain openings of all shapes and sizes, some examples of which are shown in Figure 5. At which scale do these openings begin to influence forest-atmosphere exchange? As regards light penetration, Zhu et al. (2015) propose three categories of gap sizes in temperate forests: "small" 0.49 < ∕ℎ ≤ 1 ; "medium" 1 < ∕ℎ ≤ 2 ; and "large" gaps, 2 < ∕ℎ < 3.5, where ∕ℎ is the ratio of the opening's diameter (O d ) to the mean height of the trees surrounding the gap. Openings with diameters such that ∕ℎ ≥ 3.5 are considered clearings with their own edges. Small openings, such that ∕ℎ ≤ 0.49 are not treated as gaps, because they remain in shade for much of the day. Regarding the flow of air, it is not clear at which size a canopy opening is a "pore," in that the filtering operations smooth out its effect on the flow, and at which size it is a "gap," in that it induces non-random dynamical effects. Determining a numerical threshold between pores and gaps is not possible using observations alone because of the difficulties inherent in obtaining spatially representative velocity observations in forests (Finnigan, 2000;Finnigan & Shaw, 2008), and is not straightforward even using idealized LES models, because the threshold likely to be around the scale of, or smaller than, the spatial filter. We propose, as a first approximation, that openings with horizontal diameter greater than the shear length scale (O d ≥ L s ) can be considered "gaps," because they are likely to induce fluid dynamical effects on the scale of the dominant turbulent eddies. Openings where O d ≪ L s can be considered "pores" in that they have only wake-scale effects on the flow and contribute little to the overall TKE budget (Finnigan, 2000;Raupach & Shaw, 1982;Raupach et al., 1996). The length scale L s is most soundly defined in homogeneous, dense vegetation canopies, so this is only rough approximation for patchy forests. Because L s ≈ 0.5h c , this dynamical definition of canopy gaps corresponds to the minimum size of the "small gaps" proposed by Zhu et al. (2015), which they determined with respect to light penetration rather than fluid dynamics.
Few studies have modeled the effects of patchiness on ecosystem-scale exchange. Bohrer et al. (2009) use a virtual canopy generator (Bohrer et al., 2007) to simulate three-dimensional deciduous canopies, including gaps smaller than a tree crown. The heterogeneity caused turbulent fluxes to become spatially correlated in some parts of the forests, such as stronger and more frequent ejection events occurring over shorter trees. Bohrer et al. (2009) also show that canopy gaps affect the flow differently depending on season. In winter, when deciduous forests are sparse, heterogeneity in the forest canopy causes the dominant turbulent eddies to penetrate less deeply into the canopy relative to the homogeneous case. This decreases the forest's roughness length z 0 , a common measure of surface roughness, and the displacement height h d , which reflects the height of the bulk sink of momentum. By the end of spring, however, when the PAI is higher, gaps in the forest canopy cause the dominant turbulent eddies to penetrate further into the canopy, and the values of z 0 and h d to increase. In other words, in winter, the flow perceives a patchy forest as smoother than a homogeneous one but, after spring leaf-out, it perceives a patchy forest as rougher. This probably results from there being a smaller density contrast between the gaps (high eddy penetration) and full canopy cover (low eddy penetration) in winter, as well as lower wind speeds in spring. It is not clear why Bohrer et al. (2009) focused on spring rather than summer as their leaf-out season, since, for deciduous forests, spring is a time of rapid change in terms of weather and canopy structure. The changes to the values of z 0 and h d induced by the transitions in canopy morphology are relevant to large-scale atmospheric mod-  (Hart et al., 2020). The FACE infrastructure is sited in natural gaps in the forest canopy. Note the variation in foliage color and density in (b), which was caused by insect herbivory during an outbreak of the European winter moth; (c) canopy openings in the understory layer of a tropical rainforest in Suriname; (d) canopy openings in the understory layer of a boreal forest in British Columbia, Canada. els, which employ these parameters to represent forests' effect on the atmosphere. Bohrer et al. (2009) propose that variables such as the maximum PAI and the fractional area of gaps may be used in regional models to adjust the parameterized values of z 0 , h d and the eddy penetration depth, each of which may vary by around 25% in patchy forests relative to homogeneous forests of the same density. Maurer et al. (2015) find that varying z 0 and h d seasonally, as a function of the canopy structure, produces more precise and less biased estimates of u * than models taking constant values of those parameters. In Section 6.4, we discuss the dispersive fluxes arising from canopy gaps and patchiness.

Scalar Quantities in Realistic Landscapes
Scalar processes remain far less well understood than velocity adjustment and momentum transport, in terms of both the fundamental flow statistics (Shralman & Siggia, 2000) and in geophysical applications around forests Katul et al., 2013). However, the desire to understand forest-atmosphere scalar exchange is a major motivation for measuring and modeling forests in the first place. Important scalar quantities in forest ecology include energy; trace gases such as CO 2 , water vapor, ozone, and volatile organic compounds (VOCs); ultrafine particles (UFP); litter fragments; soil particles; insect and animal detritus; and spores and pollen. In this section we focus on the numerical modeling of species that can be approximated as being passive and massless at the ecosystem scale. These include CO 2 , as well as other gases and UFPs whose lifetimes are longer than the longest air parcel residence times (Bannister et al., 2021;Janhäll, 2015;Kanani-Sühring & Raasch, 2015;Petroff et al., 2008). Other scalars, such as litter, animal detritus, and certain pollens, are much larger and must be treated differently (Section 6 below). For a scalar quantity ϕ, the conservation equation (neglecting molecular diffusion) is where is the SGS scalar flux (this term is not present when the equation is not spatially filtered, see Section 2.2 above). The term is a source/sink accounting for the emission, deposition, and production or destruction of the scalar quantity. Equation 11 can be double filtered into mean, turbulent and dispersive components, as for Equation 2. We retain the forest edge terminology from Section 3 ( Figure 4) in this section.

Scalar Transport at the Upstream Edges of Forests
Scalar quantities accumulate in the adjustment region at the upstream edge of the forest (region [ii] in Figure 4) over a streamwise distance x ≈ 9−12L c ≈ 2x A ≈ 16−20h c . Here, concentrations of the scalar quantity can be several times higher than in the air upstream of the forest. This pattern appears consistent across field observations of heat transport around forest edges (Klaassen et al., 2002), and in RANS (Sogachev et al., 2008) and LES (Kanani-Sühring & Raasch, 2015; Ma et al., 2020) models of flow around idealized forests. The turbulent fluxes of scalar quantities above the adjustment region are 1.2-3.8 times larger than the surface source rate, and are therefore compensated by horizontal fluxes from elsewhere in the forest (Kanani-Sühring & Raasch, 2015). The location of the peak scalar concentrations and turbulent fluxes is dictated by (a) the adjustment of the flow and (b) the streamwise turbulent scalar transport. The mean and turbulent fluxes are of the same order, with the turbulent component more influential in sparser forests. Strong turbulent fluxes occur above the end of the adjustment region at x ≈ x A , and z ≈ 1.5h c (Kanani-Sühring & Raasch, 2015;Ma et al., 2020). Ma et al. (2020) attribute these fluxes to scalar-rich air being swept upwards by the mean and turbulent components of the flow ( 2 is high in this region). Horizontal and vertical advection, which eddy covariance measurements do not usually account for (Aubinet et al., 2010), dominate CO 2 transport in the adjustment region (Ma et al., 2020).
Concentration peaks are more pronounced in forests with a higher mean density and are located closer to the upstream edge, for example, at x ≈ 9−12L c ≈ 8h c for PAI ≈ h c /C d L c = 8 but x ≈ 5−6L c ≈ 14h c for PAI = 2. Forest succession and management, both of which can change PAI will, therefore, greatly influence scalar concentration gradients close to edges. This is because denser forests impose more drag on the flow than sparser ones, which suppresses turbulent mixing and causes the flow to adjust more quickly, that is, the canopy drag length scale L c is reduced. Turbulent scalar fluxes are greater in magnitude and located closer to the edge with increasing forest density. The magnitude of the fluxes increases with increasing wind speed but, because the turbulent components w′ and ϕ′ vary inversely with increasing wind speed, the locations of the fluxes do not change. Higher concentrations and scalar fluxes occur at the upstream edges of forests in which the foliage is distributed uniformly than those in which the foliage is concentrated in the upper canopy (Ma et al., 2020). Forest management induces a step-change in the aerodynamics of the forest by substantially decreasing forest density. The subsequent flow pattern must then continuously adjust as the forest density increases steadily over the timescale of the management period (O(10 years)).
As to the vertical distribution of sources and sinks, scalar concentrations in the lower two-thirds of the canopy are generally higher for scalars with ground sources than for those with canopy sources (Edburg et al., 2012). The turbulent fluxes of scalar quantities are high throughout the depth of the canopy for ground sources, but their magnitudes decrease sharply toward the ground for scalars with canopy sources. Scalar quantities with canopy sources are stirred by the large turbulent eddies near the tops of the trees, whereas ground sources of scalars require intermittent turbulent motions to permeate the entire canopy depth. Scalars emitted from the canopy are therefore more evenly mixed vertically throughout the forest and have shorter residence times compared with those emitted from the ground (Edburg et al., 2012). With decreasing source height, there are larger variations in the concentrations and turbulent fluxes of scalar at the upstream edges of forests (Ma et al., 2020;. This is primarily due to low wind speeds and thus less turbulent mixing within the lower canopy layers, resulting in local accumulation of scalar concentration and large gradients. The lowest scalar concentrations coincide with the vertical locations of the sinks, if one exists (e.g., for CO 2 , in the upper canopy during the daytime in the growing season). The depth of the region of low concentration increases with distance from the upstream edge.
Around forest edges, the TKE and momentum transfer are much higher in unstable conditions than when the ABL is approximately neutral, and the flow takes longer to adjust upon meeting the canopy (Ma & Liu, 2019). In the adjustment region, the skewness of the streamwise Sk u and vertical Sk w velocity are smaller in magnitude in unstable conditions than in neutral conditions, indicating that sweep motions do not penetrate as easily into the forest canopy when the air is buoyant. In neutral conditions, CO 2 accumulates in the adjustment region at the upstream edge of the forest, and water vapor in the canopy flow region, but these patterns largely vanish in unstable conditions (Ma & Liu, 2019).
Scalar fluxes appear to adjust more slowly than momentum as the flow meets the forest, but the patterns in scalar transport are not always intuitive. Scalars are represented in models through the source term S in Equation 11 as either a concentration or flux boundary condition. In numerical models specifying flux boundary conditions, such as those discussed in the previous few paragraphs, one would expect scalar quantities to be slave to the flow because no additional length scales are introduced in the model. Based on Sogachev et al. (2008) and Belcher et al. (2012) suggest scalar concentrations and fluxes ought to reach equilibrium after a distance of x ≈ 2x A downstream of the edge, where the value of x A decreases with increasing forest density. Following this reasoning, we expect the fluxes of scalars to adjust more rapidly in denser forests. However, Kanani-Sühring and Raasch (2015) find the opposite: scalar fluxes adjust more slowly with increasing PAI, x/x A ≈ 2 for sparse forests, where PAI ≈ h c /C d L c = 1-2, and x/x A ≈ 3-4 for forests with PAI = 4.5-8. The relative size of sparse and dense patches of heterogeneity in forests may also affect scalar concentrations. For example, for flow through idealized porous media, maximum scalar concentrations occur in sparse patches where the patch size is less than the adjustment distance x A , but in the dense patches where the patch size is greater than x A (Bannister et al., 2021). The patterns in scalar transport are likely to be ecologically significant. For example, the simulations of Ma et al. (2020) indicate the air can be a few degrees warmer at the upstream edge of the forest. If this temperature variation is correct, a forest edge in the prevailing wind direction may provide very different habitats to forest only tens of meters away (Zellweger et al., 2020). Turbulent fluxes of scalar quantities adjust slowly around patchy forests and in other complex terrain. For example, turbulent fluxes of water vapor and CO 2 respectively adjust over distances x ≈ 30h c ≈ 38L c and x ≈ 60h c ≈ 76L c downstream of a forest edge (Ma et al., 2020). This amounts to some 500-1,000 m for mature forests which, given the abundance of edges in contemporary forested landscapes, indicates scalar concentrations and fluxes are out of equilibrium for much of the world's forested area.

Scalar Transport in the Lee of Forests
The area behind the leeward edges of forests shares some obvious similarities with that in the lee of windbreaks in that they are both sheltered from the flow of air. Because of the long history of windbreaks, their microme-teorology is quite established, and gives some idea of what to expect in the lee of forests, for which little observational evidence exists. The turbulent transport of scalar quantities is generally enhanced in the wake region behind windbreaks but is suppressed in the quiet zone (McNaughton, 1989), that is, for a uniform ground source, we expect higher scalar concentrations in the quiet zone and lower concentrations in the wake region, compared to the landscape mean. LES studies of scalar transport in forest lees supports this assessment. Kanani-Sühring and Raasch (2017) show that scalars accumulate in the exit region over a distance x/h c ≈ 2 downstream of the trailing edge of the forest, with peak concentrations up to around twice as high as those at a reference point far downstream. Ma et al. (2020) obtained similar results using a coupled radiation-LES model, although the magnitudes of the concentration and fluxes depended on the specification of the sources and sinks. Turbulent transport accounts for around half the total scalar transport in the exit and wake regions. The locations of the highest concentrations and turbulent fluxes do not vary with wind speed. Higher concentrations in the exit region occur at lower wind speeds, but the magnitudes of the turbulent fluxes are not affected (Kanani-Sühring & Raasch, 2017). In the exit region, the magnitudes of both the concentrations and turbulent fluxes increase non-linearly with forest density. For example, the magnitude of the turbulent fluxes for forests of PAI ≈ h c /C d L c = 8 are around two and half times those of forests with PAI = 2.

Topographical Effects on Forest-Atmosphere Exchange
Topographical effects on forest-atmosphere exchange are not random and can introduce notable horizontal fluxes that are not captured by eddy covariance or smoothed out by time averaging. In heavily populated regions, forests are overwhelmingly confined to land of marginal agricultural value, which often means sloping land (Sabatini et al., 2018). The effect of topography on forest-atmosphere interactions is therefore of wide applicability. Finnigan et al. (2020) review boundary-layer flow in hilly terrain, including sections discussing flows around vegetated hills. We refer readers to their paper for a thorough overview of topographical effects on land-atmosphere exchange. Before moving on, however, we highlight notable results from recent studies that relate to ecosystem-scale forest-atmosphere exchange of scalar quantities.

Boundary-Layer Flow Over Hills and Forest-Atmosphere Exchange
A large proportion of the literature investigating land-atmosphere exchange over complex terrain studies the flow over simplified topography, such as two-dimensional Gaussian hills or sinusoidal ridges. These studies reveal phenomena around forested hills that appear remarkably robust, despite the studies' differing treatments of the physics (B. Finnigan et al., 2020;Katul et al., 2006;. In the lee of a forested hill, an in-canopy recirculation region (ICR) may form within the forest, particularly when (a) the hills are steep, so there is a large pressure gradient at the upstream slope; and (b) in tall, dense forests, in which there is little mixing of high momentum air from above the forest with that near the forest floor. Another ICR may form at the foot of the slope at the upstream side of the hill. The ICRs result from a balancing act between the aerodynamic drag, the pressure perturbation induced by the hill, and the shear stress induced by the forests. A Figure 6. Schematic of the flow over a forested hill. A region of high scalar flux occurs around the in-canopy recirculation region in the lee of the hill. This region acts as a chimney for air parcels leaving the forest air space. helpful scale to interpret flow in hilly terrain is the hill half length, L, defined in Finnigan and Belcher (2004) as a quarter of the wavelength of the topography. The ICR upstream of the hill occurs at x/L ≈ −1 to 0, and another in the lee of the hill at x/L ≈ 2-4 ( Figure 6; B. . The air in the wake of the lee ICR is often highly turbulent (Finnigan et al., 2020). Forests absorb lots of momentum, meaning ICRs are more common in forested landscapes than in those with low vegetation, and may form in the lee of even gentle hills (Finnigan & Belcher, 2004;Patton & Katul, 2009). Both the forest density, through its effect on the canopy drag length scale L c , and the absolute height of the canopy influence the likelihood of ICRs forming. In "shallow" canopies, where ∕ < 2 2 ≈ 0.2 (Poggi et al., 2008), not all of the momentum is absorbed by the foliage and ICRs are less likely to form than for "deep" canopies, where h c /L c ≫ 0.2. Figure 6 shows the region of strong vertical fluxes of scalar quantities forms between the crest of the hill and the upstream edge of the lee ICR. In this region, scalar fluxes may be up to nearly an order of magnitude larger than the mean landscape flux, depending on the emission location and sampling height (B. Chen et al., 2020;. Scalar fluxes around the lee ICR are stronger for ground sources of scalars than for canopy sources, because the topography affects the flow and turbulent mixing more near the forest floor . In models for which multiple sources and sinks are specified, the net effect of the sources depends on the balance of the individual transport terms. B.  propose two pathways by which air parcels leave the forest: (a) the "local" pathway, where ejection events transport air parcels out the forest, approximately vertically and (b) the "advection" pathway, where parcels are transported horizontally until they encounter and are entrained into the lee ICR, before turbulent fluxes eject them from the forest. The local and advection pathways are, respectively, the dominant passages for air parcels moving from the upper and lower parts of the forest air space. The whole forest air space contributes air parcels that leave via the advection pathway, although sources from the forest floor contribute a greater proportion of the total escape. The region of high flux acts as a chimney for air parcels leaving the forest. This behavior is amenable to observational testing but has not been verified by field measurements. Zeri et al. (2010) observed higher CO 2 concentrations around a forested hill when the wind blew from certain directions, which they attributed to CO 2 accumulating in the hill's lee. However, the accumulation region fell just outside of the observational area, so the authors could not investigate this behavior further. LES simulations over realistic topography show a region of enhanced fluxes develops downstream of hill crests, near the separation point, consistent with the chimney effect that has been observed for simplified two-dimensional terrain (B. Chen et al., 2020).

Looking Forward
It is not straightforward to interpret flux measurements in complex terrain, even when armed with theory and conceptual models. Measurement towers are typically sited on the tops of hills. When considering simplified topography, this decision is not controversial. The fluxes at the crest of the hill are similar to the landscape flux, although the terrain may induce biases that cause the flux measurements to be higher or lower than the landscape mean (B. Chen et al., 2020). However, over realistic topography, while the crest usually remains the best observational location, measurements consistently underestimate the landscape flux and there is a large variability in flux measurements at different hill locations (B. Chen et al., 2020). Further testing at other sites will reveal whether this tendency to underestimate flux measurements is a common problem. Chen et al. acknowledge that they performed their experiments in neutral conditions; mixing promoted by daytime buoyancy may alleviate some of the discrepancy between the crest measurements and the mean landscape flux.
Flows induced by gravity-as opposed to pressure perturbations in the flow, as in Figure 6-are an important topic of ongoing research in forest-atmosphere exchange. In mountainous regions, gravity flows often interact with one another across time and space scales larger than the focus of this review, so we do not discuss these flows in detail here (see Oke [1988] for an introduction to gravity flows, Mahrt [1982] and Belcher et al. [2012] for related scaling analysis, and Finnigan et al. [2020] for a review of recent studies and wider discussion). Downhill (katabatic) flows can occur at night when valley surfaces emit long-wave radiation, causing the air near the surface to cool and accelerate downslope under the influence of gravity. Katabatic flows tend to decouple soil CO 2 fluxes from eddy-covariance measurements above the canopy (Aubinet et al., 2003;Feigenwinter et al., 2008), which can introduce significant errors into diurnal carbon budgets (Van Gorsel et al., 2011). In stable conditions, even gentle slopes can generate strong gravity flows (Belcher et al., 2012), meaning that, at certain sites, advection complicates the interpretation of eddy-covariance measurements at night but not during the day (Leuning et al., 2008). However, nighttime advection it is not easily addressed even with careful measurements within the forest (Aubinet et al., 2010). An interesting line of future research would be to compare low-wind nighttime eddy-covariance observations around forested hills to high-resolution LES simulations. This type of work may reveal much-needed scaling laws to determine the onset of gravity flows, the strength and even direction of which are sensitive to the balance of the inertial forces and the buoyancy (Belcher et al., 2012;Finnigan et al., 2020).

Air Parcel Residence Times
Once a molecule or particle enters a forest from outside, or is released from a leaf or the soil, it is mixed within the forest air space. This is easiest to visualize in terms of the stretching and dissipation of small air parcels. During the passage of these air parcels, the molecules or particles contained within may react, diffuse across leaf boundary layers, or deposit on surfaces. Air parcel residence times therefore affect forest ecology-for example, influencing chemical signaling (Szendrei & Rodriguez-Saona, 2010) and VOC chemical processes (e.g., Batista et al., 2019;Pugh et al., 2011), or varying the likelihood of nutrients (Fowler et al., 2009) or fungal spores (Norros et al., 2014; being deposited. The lifetimes of some VOCs emitted by forests are in the order of tens of minutes (Wolfe et al., 2011), similar to the residence times parcels moving from close to the ground. Reactive scalar quantities emitted from the ground are therefore more likely to be chemically transformed within the forest than scalars emitted in the tree crowns.
The initial vertical position of an air parcel affects its residence time, with parcels "released" near the ground having much longer residence times than those released higher in the forest canopy (Edburg et al., 2012;Fuentes et al., 2007). (Numerical models allow the flow to be studied from an Eulerian or Lagrangian point of view. The Eulerian specification of the flow focuses on specific locations in space through which the air flows with time, whereas the Lagrangian specification follows an individual air parcel as it moves through space and time. Studies looking at particle residence times in plant canopies often adopt a Lagrangian point of view.) Estimates of parcel residence times vary from a few seconds for parcels moving from the crown, to around 30 min for parcels moving from the forest floor, with spatially averaged values in the order of a few minutes (Edburg et al., 2012;Fuentes et al., 2007;Gerken et al., 2017;Hart et al., 2020). The density and morphology of vegetation influences air parcel residence times. For example, in a short, trellis-trained crop, parcel residence time increases with canopy PAI, other than for parcels released high in the canopy (Bailey et al., 2014). Residence times increase with increasing PAI, because mixing-layer-type eddies and TKE do not penetrate as deeply into the canopy . Air parcels remain in the canopy for longer when they are released from approximately the height at which most of the plant area is located. For example, for parcels released in the upper canopy, residence times are longer for top-heavy PAI profiles, such as pine plantations, than for forests with more plant area lower down. We expect longer residence times in stable atmospheric conditions, such as at night or on overcast days, because turbulent mixing is suppressed. We are not aware of any investigations into stability effects on parcel residence times in forests. Observations and RANS simulations in urban areas indicate parcel residence times generally increase with increasing atmospheric stability, although they also heavily influenced by the wind velocity and the geometry of local obstacles (Mavroidis et al., 2012).
The next challenge in this line of enquiry is to explain quantitatively how forest canopy turbulence affects air parcel residence times. As rough approximation in homogeneous forests, Gerken et al. (2017) propose that the parcel residence time (τ) is proportional to the reciprocal of the friction velocity, that is, τ ∝ 1/u * . Unfortunately, this relationship is unlikely to hold even approximately in inhomogeneous landscapes. Around forest edges, and in other complex terrain, the flow is spatially variable (e.g., Figures 3 and 6). This variation means that we cannot infer residence times across the whole forest from single variables such as u * , that is, we cannot assume a moving equilibrium (Yaglom, 1979). While there have been no investigations of the effect of forest heterogeneity on air parcel residence times, in vineyard-type canopies, heterogeneity generally decreases residence times (Bailey et al., 2014). Lagrangian investigations in urban settings show that heterogeneity induces spatial variations in the residence times of scalar quantities (Lo & Ngan, 2015, 2017 and this method may be applied to forested environments.

Part II: Improving ecosystem-scale models of forest-atmosphere exchange
In Part II we focus on recent attempts to improve numerical modeling techniques from developments in theory and experiments. Many of these topics, particularly those in Section 6, are important aspects of atmospheric science in their own right, which deserve their own reviews. However, certain techniques can be readily adapted to improve understanding of how forests and our atmosphere interact. We focus on the application of these topics in turbulence resolving models, rather than attempting to review each area exhaustively.

Drag and Plant Reconfiguration
The distributed-drag parametrization was introduced in the 1970s and remains the starting point for ecosystem-scale numerical investigations of forest-atmosphere exchange (see Table A1), usually through the approximation = − ( )| |⟨ ⟩ (Equation 7). The distributed-drag parametrization accounts for the average aerodynamic drag the plants impose over some spatial scale larger than individual twigs and leaves. For some applications, this averaging process is sufficient. Using this method, RANS and LES models can accurately resolve the mean flow through homogeneous forests and around bluff bodies (Yang, Raupach, et al., 2006;. However, the method poorly reproduces higher-order flow statistics around forests and in other vegetation canopies (Dupont & Brunet, 2008a;Ma & Liu, 2019;. The reason for this discrepancy may be lie in the fact that there is little in the formulation to account for the detailed morphology and mechanics of plants and trees. In this section, we look at numerically efficient ways to improve the drag parametrization to represent a more realistic range of forest structures. Trees and other plants reconfigure elastically to reduce drag over short time scales, visible in the everyday observation of leaves curling and tree branches thrashing in high wind (e.g., Figure 7). The reconfiguration of large trees is less extreme than of many smaller plants, which can bend elastically throughout their whole body in strong wind, significantly reducing their mean height. However, even large hardwood trees are softer geometrically and mechanically than most human-made structures, and should not be approximated as being rigid. Measurements of isolated trees in wind tunnels show that streamlining reduces the trees' frontal area by up to 54%, depending on tree species and wind speed (Rudnicki et al., 2004;Vollsinger et al., 2005). Measuring the effect of reconfiguration is more difficult in a forest than on isolated trees. Trees deep inside forests are exposed to lower wind speeds than those around edges, but they are also usually more supple than edge trees (Gardiner et al., 2016). High-frequency velocity measurements in conifer forests show that trees reconfigure at the stand scale, especially around edges, large clearings, and in the region of high shear around the top of the canopy (Dellwik et al., 2014;Queck et al., 2012). The reconfiguration of dense forests and other plant canopies is also evident in coherent "honami" waves passing through a field of wheat (Inoue, 1955;Maitani, 1979), from wavelet analysis of pine plantations (Schindler et al., 2012), and in video footage of forests taken above the crowns (Harper, BIFoR FACE site, private communication 2021).
At low wind speeds, the movement of individual leaves dominates foliage reconfiguration (Tadrist et al., 2018), for example, curling up to form cones and cylinders (Vogel, 1989). This foliage-dominated configuration could account for the scatter observed in the measured values of coefficient C d at low wind speeds. At low speeds, the total plant area exposed to the wind varies more than at higher speeds; at low speeds, some leaves will be sheltered only until a gust rolls though the forest. Measuring drag coefficients using instantaneous velocity measurements-that is, calculating the drag as being proportional to | | rather than U 2 -reduces scatter in estimations of C d and the wind-speed dependency of the product C d a(z) (Cescatti & Marcolla, 2004;Queck et al., 2012).
Plant reconfiguration is particularly relevant in patchy forests because the flow is constantly adjusting and readjusting as it moves across gaps and clearings. Capturing this reconfiguration directly in ecosystem-scale models is currently unworkable because of the computational expense of simulating the motion of flexible bodies in turbulent flow. However, it may be partly captured by modifying the drag parametrization in Equation 7. For example, one could consider the PAD, represented in a(z) in Equation 7, as varying with velocity (Speck, 2003), or replace the drag coefficient C d , with a shape factor that varies with the flow (Gaylord & Denny, 1997). However, these approaches require a detailed three-dimensional knowledge of the forest structure and its response to wind loading. A more empirical approach is to proceed from the observation that reconfiguration leads to a lower increase of drag with velocity than could be expected by assuming a(z) and C d both take constant values (Vogel, 1989). Therefore, instead of having the drag ∝ 2 , we have ∝ 2+ , where B is known as the Vogel number (De Langre, 2008;Vogel, 2020). B therefore acts as an empirical modification to Equation 7 to account for the variation of a(z) and C d with velocity. Negative values of B imply that, because of reconfiguration, plant drag increases with velocity more slowly than the usual assumption of velocity squared. Where the plants do not reconfigure, B→0 and the drag increases approximately with the square of the velocity. Pan, Chamecki, and Isard (2014) use LES to model the reconfiguration of a maize canopy. They maintain the quadratic dependence of drag on velocity in the drag parametrization in Equation 7, accounting for plant reconfiguration using a velocity-dependent drag coefficient = (| |∕ ) , where A v is a velocity scale related to the shape and rigidity of the plants. Using the variable C d , their LES model simulated flow with more accurate higher order statistics than the same model with a constant C d value (with similar performance with respect to the mean velocities, momentum transfer and TKE). Using a variable C d value also improved estimates of momentum transport by sweep motions penetrating the canopy. In similar work, Pan, Follett, et al. (2014) show plants reconfigure more at higher flow velocities, reflecting the results of lab investigations on isolated plants (Tadrist et al., 2018). In both Pan, Chamecki, and Isard (2014) and Pan, Follett, et al. (2014), the absolute values of Sk u and Sk w increase with a more negative Vogel number-that is, simulating more reconfiguration-as did the ratio of sweeps to ejections. To see why, taking B = −0.74 from Pan, Chamecki, and Isard (2014) gives ∝ 1.26 in Equation 7 rather than ∝ 2 . The drag force therefore increases more slowly with increasing velocity, leading to a greater contrast of TKE inside and above the canopy. The decrease in the velocity exponent also increases the frequency and strength of sweep motions because strong events (u ' > 0) penetrate further into the canopy. The canopy drag length scale L c increases when the value of C d is smaller because stronger events are able to penetrate further into the forest before drag halts their progress.
Field studies of plant reconfiguration typically focus on time-averaged adaptation over long time intervals, although measurements in poplar crowns show C d values decrease with increasing wind speed across averaging periods as short as 1 s (Koizumi et al., 2010). Measurements in a patchy conifer canopy found B ≈ −0.8, although this value varied with the wind direction . Pan, Follett, et al. (2014) ran simulations using B = 0, −2/3, −1, and −4/3, finding that a non-zero Vogel number improved the agreement of the simulations with flume observations, with the greatest fidelity obtained using B = −1. Other estimates for B include −1 < B < −2/3 from biomechanics theory . B ≈ −1 appears to be a robust measurement for poroelastic structures, including forests and other vegetation canopies (Gosselin, 2019 and references therein). We suggest 0 and −1 as approximate upper and lower bounds for the Vogel number B.

Terrestrial Laser Scanning
Most studies assume that the local foliage density a(z) varies only with height, if at all (Table A1). For example, Figure 8a shows a(z) derived using the parametrization by Lalic and Mihailovic (2004), as employed by Yan et al. (2020), among many others. The forest is assumed to be horizontally homogeneous at each height, reducing the forest morphology to a single dimension (z), typically with most of the density confined to the tree crowns. This is a reasonable approximation for certain forests, such as unthinned conifer plantations (Figure 8c). However, it is a poor approximation for many forests, such as oak-dominated temperate forests (Figures 5 and 7b) and certain tropical forests (Figure 9), which have dense understory layers in which the PAD is similar to that of the crown layer, or multiple modes of leaf density (Schneider et al., 2019;Zhao et al., 2011). Most models also assume that forest edges are structural continuations of the main body of the canopy (Dupont & Brunet, 2008a;Kanani-Sühring & Raasch, 2015;Ma et al., 2020). Unless the forest has been recently disturbed, this assumption is unrealistic because many trees and understory plants grow into the light to maximize the leaf area available for photosynthesis, thereby forming an almost closed edge (Figure 8b).
A major technological leap since the first LES investigation of forest-atmosphere interactions (Shaw & Schumann, 1992) has been the development of terrestrial laser scanning (TLS), also known as terrestrial LiDAR. TLS provides non-destructive, in situ information of the three-dimensional structure of trees and forests, down to a resolution of millimeters (Calders et al., 2020). See, for example, a captivating example of a TLS scan of an Australian notophyll vine forest in Figure 9 (Calders et al., 2020). Note the variety of plant forms, stem diameters, and leaf sizes for this uneven-aged, biodiverse, forest. Note also, in the right-hand panel of Figure 9, the pronounced increase in PAD at about 2 m, suggestive of a herbivore "browse line" (Hazeldine & Kirkpatrick, 2015), the rel-  atively sparse canopy between 2 and 10 m, and the slow drop-off of PAD with height above about 28 m, typical of a biodiverse forest. TLS is currently used for a range of ecological applications (Ashcroft et al., 2014;Momo Takoudjou et al., 2018), but has yet to be adopted widely in models of forest-atmosphere exchange.
TLS combined with LES has the potential to revolutionize the realism of ecosystem-scale models of forest-atmosphere exchange. TLS allows real sites to be reconstructed numerically by splitting the forest volume into cells, each of which is assigned a plant density value calculated from TLS point clouds (Boudreault et al., 2015;Calders et al., 2015;Queck et al., 2012;Raumonen et al., 2015). The few studies that have used this approach in heterogeneous forests show important flow behavior may be missed by assuming heterogeneity. Schlegel et al. (2012Schlegel et al. ( , 2015 use TLS measurements in an LES model of flow around a large forest gap. They show that small-scale plant heterogeneity creates sustained upwards motion in denser patches of forest and downwards motion in clearings and large gaps. Failing to account for the true structure of forest edges strongly overestimates the penetration of air into the forest (Schlegel et al., 2012(Schlegel et al., , 2015, and neglects the possibility of dispersive fluxes of momentum and scalar quantities (Boudreault et al., 2017;Li & Bou-Zeid, 2019). TKE and Reynolds stresses decay faster behind closed edges and strong cross flows may develop. Away from the edges, large gaps and clearings deflect the flow downwards, creating advective fluxes within the forest air space (Queck et al., 2016). Boudreault et al. (2017) use TLS measurements (outlined in Boudreault et al. [2015]) and LES to model air flow across a forest edge. Compared to a homogeneous edge, the gaps induce variations in the flow, for example, the streamwise (U) and vertical (W) velocity components respectively vary by around 20% and 5% of their spatial means at z = 0.5h c . The turbulent momentum fluxes are likewise higher in the heterogeneous case, for example, σ w varies by around 40% of its spatial mean, because air is forced through patches of low density.

Stochastic Drag Forcing
The product C d a(z) in Equation 7 is typically approximated as a smoothly varying function of height, with C d taken as constant. However, Finnigan and Shaw (2008) raise the important point that the dominant large eddies around forests have diameters in the order of L s ∝ h c (Bailey & Stoll, 2016;Raupach et al., 1996). To resolve the structure of these eddies in numerical models, the vertical filter needs to be much smaller than h c , for example, a filter of h c /25 is taken in the model in Appendix A. At this resolution, we can no longer assume a(z) is a smooth function-for example, notice the structural variation of the forest understories in Figures 5 and 8. Finnigan and Shaw (2008) propose representing C d a(z) through a stochastic variation overlaying a smooth background trend, thereby introducing a stochastic forcing into the resolved momentum equations.
To illustrate this argument, we consider a small, horizontally homogeneous forest patch, represented using the drag parametrization in Equation 7, with the transport equations solved using LES (see Appendix A for numerical details). We define two cases, with the same dimensions and mean density. For Case 1, C d a(z) is uniform throughout the domain. For Case 2, C d a(z) varies randomly in space throughout the forest by ±2% of the value for Case 1, a small variation in comparison to the natural variation of forests. The random spatial variation of C d a(z) is almost imperceptible in the mean variables and second-order statistics such as the TKE and kinematic turbulent momentum flux. However, Figure 10a shows the streamwise velocity skewness Sk u at the upstream edge of the forest (x/h c ≈ 0-5) is less negative for Case 2 than it is for Case 1. The value of Sk u above the forest is more negative for Case 2 than it is for Case 1. The streamwise velocity kurtosis Kt u at the upstream edge of the forest (x/h c ≈ 0-5) is smaller for Case 2 than it is for Case 1 within the forest, and larger for Case 2 than it is for Case 1 above the forest (Figure 10b). These statistics indicate the small stochastic forcing in Case 2 decreases the coherence of the flow at the upstream edge-that is, the upstream edge behaves less like a bluff body than the homogeneous Case 1-leading to fewer lulls in the streamwise velocity within the forest but more lulls just above the canopy. We do not generalize these results further because the cases are idealized. But they support the conclusions of Dupont and Brunet (2008a) and Pan, Chamecki, and Isard (2014), for example, that the probability distribution of gusts indicated by higher-order statistics are sensitive to the model setup, including any stochastic element. This sensitivity needs careful testing.

Waving Plants and Biological Backscatter
The velocity spectra of airflow around forests and other vegetation canopies often contain peaks that correspond to plant movement (Cava & Katul, 2008;Dupont et al., 2018;Finnigan, 1979aFinnigan, , 1979b. These peaks indicate energy is transferred in two directions-the flow perturbs the plants, but the plants also perturb the flow. The latter effect is usually ignored in models at the ecosystem scale and above on the assumption (usually implicit) that the turbulent structures generated by plant movement are much smaller and less energetic than the dominant mixing-layer eddies. However, this neglects the possibility of resonant effects occurring when the passage frequency of the dominant eddies approaches the natural frequency f 0 of the moving plants, as has been observed in crops (Py et al., 2006). Trees near the edges of forests, such as in patchy landscapes and around clearings, are more susceptible to resonance effects than those further inside forest stands (Dupont et al., 2018).
Accounting for this two-way transfer of energy is not straightforward at the ecosystem scale. Dupont et al. (2010) use LES to model a crop canopy as a poroelastic continuum, with plant movement incorporated into the momentum equations as small mechanical oscillations of rigid stems. Other researchers had previously developed models coupling wind flow and plant swaying in a similar way, but with analytical solutions obtained from simplified velocity statistics rather than through LES (Gosselin & de Langre, 2009;Webb & Rudnicki, 2009). Pivato et al. (2014)  The direct approaches of Dupont et al. (2010) and Pivato et al. (2014) require the plant architectures to be heavily simplified to be computationally feasible. This is not a major problem for tall, slender trees such as maritime pine, the subject of Pivato et al. (2014). However, decurrent trees, which include many broadleaf species, are structurally more complex and can have modes of vibration across several scales, for example, f 0 ≈ 0.5 Hz in the trunks and several Hz in the branches (Schindler et al., 2013). Capturing these interactions at the scale of an entire forest would be very demanding computationally. A possible workaround in patchy landscapes is to proceed from the observation that wind velocity spectra at forest edges contain peaks around f 0 , the natural frequency of the trees (Dupont et al., 2018). This motion is especially visible in the spanwise direction because the turbulent velocity perturbations are smaller than those in the streamwise direction. From a modeling point of view, the trees' movement transfers energy from the SGS to the smallest resolved scales (Piomelli et al., 1991), a process known as "backscatter." In general, backscatter is most apparent where small but energetic eddies are present (Mason & Thomson, 1992), such as around forests and other complex structures. Studies of urban canopy flow provide a template of how backscatter could be incorporated into forest models. For example, O'Neill et al. represent the stochastic effects of backscatter in their LES simulations of the neutral surface layer (O'Neill et al., 2015) and street canyon flow (O'Neill et al., 2016) by incorporating random acceleration fields a i in the momentum equations. This gives where ν SGS is SGS eddy viscosity and the ellipsis represents the other terms carried over from Equation 2b. Here, the Smagorinsky SGS scheme is shown. In principle, this approach could be adopted for an LES of forest canopy flows, ideally with the acceleration terms a i spaced at frequencies corresponding to the movement of tree parts, while ensuring zero divergence. Selino and Jones (2013) adopt a similar approach for a very different purpose, using synthetic SGS turbulence to improve video animations of trees moving in wind.

Resolved Trees
The difficulties involved in incorporating realistic structure into the distributed-drag parametrization, discussed above, raise the question of whether it is possible to model the forest elements directly. Poggi and coworkers approximate plants as rigid circular cylinders in water-flume experiments, finding the wakes in the lee of the vegetation stems perturbed the dominant mixing-layer type eddies at the crown top (Poggi & Katul, 2008a;Poggi et al., 2008;. Yue et al. (2007) model corn plants as stems surrounded by leaves, applying a distributed-drag parametrization at the leaf points and, at the stem points, a force calculated as drag around a cylinder.  and  use LES to investigate the flow around isolated, solid tree-like fractals, finding that model's predictions of the drag imposed by the fractals were sensitive to the orientation of the branches. Böhm et al. (2013) modeled trees as cylindrical trunks below spherical crowns in a wind tunnel. They conclude similarly to Poggi, Porporato, et al. (2004) that the partially coherent wakes of the bluff canopy elements modify the dominant eddies at the crown top, and that the flow within the canopy divides into wake and non-wake regions, where the wake regions scale with the stem diameter (Ghannam et al., 2020). Yan et al. (2017) performed high-resolution LES studies of a regular array of bluff elements, specified similarly to the wind-tunnel model of Böhm et al. (2013). In a series of separate simulations, Yan et al. (2017) configure the trees as (a) entirely bluff bodies, (b) solid trunks with the crown represented as distributed drag, and (c) as entirely distributed drag. The spatially averaged flow statistics were similar across the three cases, with slightly higher turbulence in the shear layer using the bluff-body representations. As an interesting alternative, Schröttle and Dörnbrack (2013) use LES to simulate flow around 16 Pythagoras trees (A type of fractal constructed iteratively from a right-angled triangle with squares erected on each of its sides), treated as immersed boundary layers with the outer tree branches 3.35 K warmer than their surroundings. They show the thermally driven vortices from the trees, of diameter roughly h c and turnover time of ∼30 s, interfere with the shear-generated coherent structure at the tops of the trees. However, the study was a method prototype, and its results reveal little about real forests.
For the time being, because of the computational expense of simulating turbulent flow, resolving forest structure directly remains out of reach for ecosystem-scale investigations of forest-atmosphere exchange. For example, Yan et al. (2017) use an extremely high resolution model to discretize the trees (dx = dy = 0.03h c ), which severely constrains the number of trees that can be simulated. Further, treating the trees as bluff bodies does not account for plant reconfiguration, may overestimate the scales of vortical wakes behind tree crowns, and excludes the possibility of resolving the cooperative waving motion that occurs in real vegetation. There are no theoretical constraints on the number of trees that can be included in wind-tunnel models, although practical considerations to preserve scale relationships may of course impose limits. It is also difficult to include canopy exchange and other ecological processes in physical models of forests, and to maintain flow with sufficiently high Re values, because higher values are needed than in many buff-body engineering applications (Gromke, 2018). Results from bluff-body models may be useful to derive drag parametrization schemes for use in larger scale simulations. Böhm et al. (2013) motivate this approach by identifying that wake-scale TKE in bluff-body canopies is around 1/5 L s rather than 1/100 to 1/10 L s typical of vegetation canopies. This implies momentum is transferred more efficiently to vegetation than to canopies of bluff bodies. Further, information theoretical analysis of flume experiments suggests the von Kármán streets formed in the wake of stems can cause scalar variance to be transferred from small to large scales (Ghannam et al., 2020), that is, opposing the usual turbulent cascade characterized by constant fluxes of energy and scalar variance from the scales of production down to dissipation.

Modeling In-Canopy Chemistry
BVOCs are ecologically important in forests (Niinemets, 2010;Visakorpi et al., 2018) and influence air quality, meteorology, and the climate through their interactions with oxidants such as O 3 and OH (Lelieveld et al., 2008;Peñuelas & Staudt, 2010;Rap et al., 2018). Plants release a wide variety of compounds from tissues above and below ground. Flowers and fruits release the widest variety of compounds, while leaves generally have the greatest mass emission rates (Laothawornkitkul et al., 2009). The most important class of BVOCs for air quality and climate is the terpenoid class of compounds, particularly the five-carbon (C5) hemi-terpene, isoprene, the C10 monoterpenes such as α-and β-pinene, and the C15 sesquiterpenes such as β-caryophyllene. Leaf-scale emissions of isoprene depend sensitively on temperature and light (and, hence, forest structure). Emissions of terpenoids vary exponentially with leaf temperature (Guenther et al., 2006). This class of compounds displays a great range of lifetimes with respect to reaction with OH and O 3 , from several tens of minutes for isoprene for reaction with OH (Apel et al., 2002) to 2 min for the reaction of β-caryophyllene with O 3 (Shu & Atkinson, 1994). For reaction, BVOC molecules must be well mixed at the molecular scales. Because the source of oxidants (i.e., OH and O 3 ) for the degradation of BVOCs is outside the forest, imperfect mixing acts as an apparent brake on reaction rates, a process known as segregation Krol et al., 2000;Pugh et al., 2011). Fluxes of unreacted terpenoids out of forest patches therefore depend sensitively air parcel residence times, which vary with position within the canopy, as discussed above in Section 4.4.
Studies investigating BVOC chemistry in forests are usually interested in time and space scales that are too large for the turbulent flow to be resolved by DNS, LES or RANS, requiring the turbulent exchange to be parametrized Forkel et al., 2015;B. Wang et al., 2017). These parametrizations are typically based on K-theory, which should be used cautiously in forests and other vegetation canopies (Finnigan, 2000;Monteith & Unsworth, 2008). K-theory works well in situations where the sum of the flux transport and buoyancy terms of the scalar flux budget is small compared to the production term (Bash et al., 2010;Freire et al., 2017)-for example, in neutral conditions, where there is a strong vertical gradient in the scalar concentration, such as investigating isoprene released from leaves. However, the gradient diffusion assumptions in K-theory often fail in real forests-for example, in stable atmospheric conditions or when considering scalars with multiple and temporally varying sources and sinks. Model predictions of BVOCs and their oxidation products can be very sensitive to the turbulence parameterization used (Bryan et al., 2012;Makar et al., 2017). It is therefore important to make the parametrizations as robust as possible.
Ecosystem-scale LES models can simulate counter-gradient transport and are therefore likely to be indispensable tools in developing chemistry parametrizations that scale to large space and time resolutions. Due to the complexity of the task, there have been few attempts to investigate forest chemistry while resolving turbulence. However, the urban literature is an excellent source of relevant techniques (e.g., Bright et al., 2013;Buccolieri et al., 2018;Khan et al., 2020;Kwak et al., 2015;Liao et al., 2014;Zhong et al., 2016, and references in each). One possible path is to couple chemistry models to LES, a technique which has recently been used to investigate chemical transformation, transport, and deposition of air pollutants in realistically shaped urban areas (Khan et al., 2020). However, the computational expense of the coupled approach heavily restricts the resolution of the domain and the number of chemical species that can be investigated in each run. Another approach is to model chemistry inside a forest air space using "box models," which treat the air space as a fixed volume into which species are emitted and are able to react. Box models require the characteristics of the turbulence to be specified a priori, such as through an exchange velocity between the box and its surroundings. This simplification allows computing resources to be reallocated to more complex chemistry or particle microphysics than is possible when the turbulence is highly resolved using RANS or LES. Box models are widely used to model street canyon chemistry (Holmes & Morawska, 2006;Murena, 2012;Zhong et al., 2016) and have been used for one-dimensional investigations of forest-atmosphere exchange across large homogeneous canopies Pugh et al., 2011). In patchy landscapes, multiple boxes would be required to account for the fluid dynamical regions that form around forest edges and clearings. The urban literature offers a precedent for dividing air spaces into dynamical regions. For example, models of street canyon chemistry have used multiple boxes to represent the "compartmentalization" of fluid dynamical phenomena such as counter-rotating vortices (Dai et al., 2021;Kwak & Baik, 2014;Murena, 2012;Zhong et al., 2016Zhong et al., , 2018.

Modeling Particle Deposition
Ecosystem-scale investigations of forest-atmosphere interactions are often motivated by questions concerning trace gas exchange, usually representing the species of interest as passive scalars (see Section 4). However, the behavior of many biologically important particles cannot be approximated in this way. For example, many pollens and spores have substantial mass and are subject to inertial forces different to those on a trace gas molecule (see, e.g., Hinds, 1999;Seinfeld & Pandis, 2016). Freshly nucleated UFPs, resulting from the oxidation of BVOCs (Kulmala et al., 2001(Kulmala et al., , 2007, are produced in high number concentrations around forests and may coagulate (Kulmala et al., 2001;Pierce et al., 2012). These processes introduce physics requiring a different mathematical approach to that for fluid flow or particle growth/evaporation (Jacobson, 2005;Seinfeld & Pandis, 2016;Spracklen et al., 2006).
Many studies of scalar quantities  and particles (Aylor & Flesch, 2001) in vegetation adopt a Lagrangian stochastic modeling (LSM) approach, which requires the vertical profiles of turbulence statistics to be specified a priori. However, around patchy forests, determining a priori vertical profiles is extremely difficult because the dynamics are so spatially varied. There have been a handful of attempts to incorporate non-Gaussian turbulence in an LSM framework, with varying success (e.g., Reynolds, 2012). More recently, several groups have adopted an Eulerian specification of the flow field to model particle deposition. Around patchy forests, Eulerian models offer the advantage of resolving the velocity statistics directly down to the scale of their grids. This does not necessarily affect the predictive ability of the models around vegetation. For example, Gleicher et al. (2014) simulated more accurate concentrations of spores in a homogeneous maize canopy using an Eulerian LES model than using their equivalent LSM. However, an inherent difficultly using LES is that canopy deposition occurs at spatial scales much finer than the spatial filter and therefore must be parametrized. One method is to include a sink term in the conservation equation, where α is the attenuation coefficient-which accounts for the flow's response to the forest density (Cionco, 1978)-ϕ l the resolved local particle concentration, and E ϕ the efficiency of particle deposition (Friedlander, 2000;M. Lin et al., 2012;X. Lin et al., 2018). Pan, Chamecki, and Isard (2014) use a similar approach, modifying the deposition model in Aylor and Flesch (2001) to generate a sink term linked to PAD, This formulation considers the distribution of the forest density directly though the incorporation of PAD, as a(z), and the projection coefficients P x and P y , which respectively account for the PAD facing the streamwise and spanwise directions. The ground deposition of particles can be modeled using a surface flux boundary condition. X. Lin et al. (2018) and Pan, Chamecki, and Isard (2014) use slightly different formulations for the deposition efficiency E ϕ based, respectively, on a parametrization of molecular diffusion and on observations of particle impaction onto cylinders. The deposition efficiency term E ϕ accounts for the momentum and size of the particles, so that Equations 13 and 14 are solved separately for each particle size, with the only difference between the solutions resulting from the approximation of the deposition efficiency. The same is true for the surface-flux parametrization representing ground deposition.
The mechanisms of particle deposition, and hence deposition velocities, are highly size-dependent for the size ranges of particles commonly encountered in forests (Litschike & Kuttler, 2008). X. Lin et al. (2018) found the deposition velocity decreased sharply with increasing particle size for the range of sizes they investigated (diameters 10-50 nm). Pan, Chamecki, and Isard (2014) investigate size indirectly through the ratio of the particle settling velocity w s to the friction velocity u * , with "light" particles having w s /u * ≈ 0.04 ≪ 1, and "heavy" particles w s /u * ≈ 0.2. They show, for large values of w s /u * (heavy particles) and sources near the ground, few particles escape the canopy, reflecting their empirical observation that plant diseases initiated deep within a canopy move upward to the canopy top before spreading. Unsurprisingly, more particles are deposited in denser vegetation (Branford et al., 2004), particularly at z ≈ h c , and therefore lower concentrations of particles occur deep in the canopy. The rate of dry deposition generally increases with increasing turbulence because the thickness of the quasi-laminar boundary layers around plant elements is reduced, therefore increasing the probability of impaction (Fowler et al., 2009). However, the rate at which particles are deposited depends on the tree species. Studies of urban trees indicate that particles deposit onto conifers more efficiently than onto broadleaf trees (L. Chen et al., 2017;Pace & Grote, 2020). Beyond these general observations, there appears to be no clear dependence of fluxes or deposition of fine particles on broad-brush measures of the canopy morphology (Katul et al., 2011;X. Lin et al., 2018). Deposition patterns appear to depend strongly on the complex arrangement of plant area in forests and other plant canopies (Fowler et al., 2009).
A future task, straddling in-canopy chemistry and trace-gas deposition, is to investigate how forest patchiness and more complicated weather conditions affect the rate of deposition and stomatal uptake of trace gases. For example, particle impaction can be several times higher in unstable conditions versus stable ones (Chiesa et al., 2019;Fowler et al., 2009;Pryor et al., 2008). This suggests leaf related removal by forests is likely much lower at night than measurements taken during the day would suggest; impaction is much lower in stable, nocturnal conditions and the stomata are mostly closed at night. Local atmospheric stability gradients, such as from patchy heating from sunlight, may produce apparent fluxes when particles and gases trapped in stable conditions are eventually released in convective plumes (Chiesa et al., 2019). Most models of spore dispersal in forests consider only dry deposition, which is a major simplification in many climates where forests are found, and pay little attention to the resuspension of deposited particles. However, the rates of particle removal by rain and resuspension by depend on species composition and the local meteorological conditions, such as the frequency and intensity of rainfall events (L. Chen et al., 2017). Airborne fungal spore concentrations in forests are generally higher in wet conditions (Crandall & Gilbert, 2017), which are physiologically favorable to certain fungal species such as Ganoderma spp. (Stępalska & Wołek, 2009). Air vortex rings, which can carry dry-dispersed spores away from the host plant, form around the impact site when raindrops hit plant surfaces, (Kim et al., 2019). This process, which is not accounted for in current models, is likely significant in the transmission of fungal spores because particles transported by the air vortices can reach beyond the laminar boundary layer around plant surfaces, enabling long-distance transport in the turbulent flow.

Heterogeneous Sources and Sinks of Scalar Quantities
Section 4 considered investigations where scalar quantities were approximated to have horizontally homogeneous sources and sinks. This approach provides a good conceptual template for certain processes around forests, such as the encroachment of pollutants from surrounding areas, or the release of isoprene, which is overwhelmingly from sunlit leaves in the canopy from a monoculture stand (Sharkey et al., 1996). However, in mature forests, most scalars have multiple sources and sinks, whose distribution varies temporally and spatially. This heterogeneity often requires dense networks of equipment to measure, and is therefore not widely reported, despite the known difficulties in generating spatially representative observations in forests (Aubinet, 2008;Finnigan, 2008;Massman & Lee, 2002). Visible examples of this heterogeneity include the water vapor plumes that form over forests after it has rained, for example, during periods of high convection but a low lifting condensation level height over tropical forests (Jiménez-Rodríguez et al., 2021). Less obvious examples include the steep temperature gradients that can form within forests on clear days. In direct sunlight, broad leaves may be 20°C warmer than their surroundings (Monteith & Unsworth, 2008;Schuepp, 1993;Vogel, 2009) and can therefore act as highly localized heat sources.
We are not aware of numerical investigations of how patchy sources and sinks affects scalar transport within forests. LES of unevenly heated generic surfaces suggests thermal heterogeneities drive the local mean flow in certain conditions, such as when the geostrophic wind is weak. In these conditions, the dispersive fluxes of heat account for more than 40% of the total sensible heat flux at z = 100 m and up to 10% near the surface (Margairaz et al., 2020). The effect of differential heating is not easy to constrain in models of forests because the locations of sunflecks (brief periods of high photon flux density) change quickly depending on branch movement (Way & Pearcy, 2012). The radiative properties of the leaves may even vary between the sun and shade sides of a tree (Vogel, 1968). One possible approach is to combine the modeling techniques of thermal transfer in urban areas (e.g., Martilli et al., 2002;Salamanca et al., 2010) and radiation transfer in vegetation (Ma & Liu, 2019). For example, the leaves of the forest may be represented using a probability density function (PDF) of small, flexible, surfaces with high absorptance, and the trunks through a PDF of vertically aligned cylinders of varying thickness and low absorptance. This type of testing is probably best carried out in uniform canopies to begin with, so as to differentiate between patterns in scalar transport resulting from terrain effects and those resulting from spatially heterogeneous sources and sinks. For uniform sources around forest edges, the patterns of scalar concentrations and fluxes are most complicated in the adjustment region where they are influenced by the strong turbulence, the inflow concentration from the surrounding environment, and the vertical distribution of the sources and sinks. Eventually, however, it may be possible to investigate the effect of heterogeneous source and sinks on scalar transport in patchy Section 7.1 revisits this point again in the context of customized numerical models of real study sites.

Calculations of Dispersive Fluxes
Acknowledging the spatial heterogeneity of real forests raises the question of how to handle the dispersive fluxes of momentum and scalars that result from the double averaging of the momentum and transport equations. The dispersive fluxes are the spatial correlations in the time-averaged statistics; the dispersive flux of momentum, for example, is given by the term − ⟨̄̄⟩∕ in Equation 2b. The dispersive fluxes are usually excluded from analyses of forest-atmosphere exchange (Kaimal & Finnigan, 1994;Patton & Finnigan, 2012), which is a reasonable assumption in dense, homogeneous canopies in which the dispersive fluxes of momentum are usually small (Moltchanov et al., 2011;Poggi & Katul, 2008b;. However, in sparse canopies-roughly comparable to the densities found in savannahs-dispersive fluxes can be ≈ 30% of the conventional momentum flux, even when the canopy is homogeneous (Poggi & Katul, 2008b). Further, investigations in model canopies show structural inhomogeneities can generate large dispersive fluxes of momentum (Harman et al., 2016;Moltchanov et al., 2011Moltchanov et al., , 2015. Boudreault et al. (2017) show that spatial variability in forest canopy structure induces dispersive fluxes that account for 10%-70% of the total variances of U and W at the upstream edge of the forest, across a streamwise distance of x/h c = 0-8 ≈ 0-x A . Here, the dispersive momentum fluxes and skewness are greater than their turbulent counterparts, for example, with the dispersive flux of momentum accounting for >50% of the total momentum flux. Away from the edges, gaps and other patchiness decrease the efficiency of momentum transfer at the crown top. This suggests gap-induced flow phenomena interfere with the mixing-layer-type coherent structures that form around the tops of the trees. Bailey et al. (2014) observed a similar result in row-structured versus homogeneous trellis-trained crops. This interference is probably strongest when the structural inhomogeneities are of a size ≈ L s , so that the vortices shed around them are of a similar scale to the coherent structures that develop around the crown top. In a more idealized example, Li and Bou-Zeid (2019) use LES to show that rough, heterogeneous surfaces affect the dispersive fluxes of momentum and scalar quantities more than the turbulent components of the fluxes. As in Boudreault et al. (2017) and Li and Bou-Zeid (2019) show the dispersive components can comprise most of the total momentum flux.
Wind-tunnel measurements in homogeneous plant canopies found the dispersive scalar fluxes to be small . However, dispersive fluxes of scalar quantities generated by patchiness or spatially heterogeneous source/sinks have not been studied in real forests. Recent studies in urban areas and generalized porous media suggest the dispersive fluxes of scalar quantities should not be dismissed out of hand. Philips et al. (2013) use LES in an urban canopy to show that the dispersion of a scalar quantity is sensitive to the geometry of the obstacles surrounding the source (they observe a plume's evolution directly, rather than investigating time-averaged quantities). Li and Bou-Zeid (2019) show the dispersive fluxes of scalar quantities do not always follow the flow of momentum, with obstacle geometry typically affecting dispersive fluxes of momentum more than those of scalars. The authors attribute this difference to the physical mechanisms involved. The air's velocity decreases before it reaches the upstream face of an obstacle, and therefore pressure affects momentum transfer away from surfaces. However, the air must touch an obstacle's surface to deposit or take up scalars, so scalar transport is much more spatially confined. LES studies of ABL flow over homogeneous surfaces indicate dispersive fluxes of heat are modulated by two broad flow regimes. The first is where surface heterogeneities, such as unevenly heated ground, drive the dispersive fluxes. The second is where dispersive fluxes are driven by turbulent coherent structures in high-shear conditions (Inagaki et al., 2006;Margairaz et al., 2020). Numerical investigations of patchy forests should include calculations of the dispersive fluxes. Ignoring these fluxes may mean ignoring over half of the total momentum flux, for example, by focusing only on deviations from the time average and neglecting spatial deviations. Real forests channel air into gaps and patches, creating spatially coherent structures whose contributions can be as large as those induced by turbulence.

A Note on the Resolution and Domain Size in LES Models
Scientists continually face a trade-off between scale and resolution. This choice is particularly relevant to fluid modeling because of the extreme computational expense of simulating turbulence. LES of ecosystem-scale interactions around forests currently use around 10 6 -10 7 cells (e.g., Ma et al., 2020), although Patton et al. (2016) deployed some 4 × 10 9 cells in an enormous computational effort. Researchers must decide how to deploy those cells most effectively. Two competing demands must be balanced: (a) the LES domains should be large enough to simulate the largest boundary-layer eddies; but (b) should have a filter width fine enough to resolve the smallest eddies sufficiently (Wood, 2000). As to the first requirement, we denote the domain size (L x , L y , L z ), where L i is the domain length in the ith direction, and z i the boundary-layer depth. For neutrally stratified conditions, a domain size of at least (≈3z i , ≈3z i , z i ) is probably sufficient (Mason & Thomson, 1987;Wurps et al., 2020). Larger horizontal domains-e.g. (8z i , 8z i , z i )-are needed to capture the largest eddies in free convection conditions. In stable conditions, smaller domains are required in an absolute sense, because the largest ABL eddies scale with z i , which may be 100 m or less, compared with z i ≈ 1-2 km in highly unstable conditions (Basu & Porté-Agel, 2006;Wurps et al., 2020). However, using too small a value of L z , particularly over obstacles, artificially negates any interaction of the boundary layer and the residual layer aloft (Grylls et al., 2020); L z ≥ 3z i appears to be a reasonable minimum for stable conditions (Basu & Porté-Agel, 2006;Wurps et al., 2020).
Determining an appropriate filter resolution for LES requires careful consideration of the atmospheric conditions, the forest structure, and the topography. A useful rule-of-thumb for ABL investigations is that a filter width Δf < z i /60 is sufficient for the first-and second-order flow statistics to be well resolved . This equates to Δ f ≈ 2.5, 10, and 20 m for stable, neutral, and convective conditions, respectively, over unvegetated ground (Wurps et al., 2020). However, for LES models of forests, we need to consider the dominant turbulent eddies in the shear layer at the top of the canopy. In neutral conditions, the length scale for these eddies is ≈ 2 2 ∕ ( ) (by rearranging Equation 8), giving L s ≈ 7 m for β = 0.35, C d = 0.2, and a(z) ≈ 0.18 m 2 m −3 , typical daytime values for mature forests. To resolve the eddies' structure, we need the vertical filter width to be significantly smaller than L s , say L s /8. A vertical filter width within the canopy of around a meter or less is a good first approximation. The dominant eddies have a mean streamwise spacing of around 8.1L s ≈ 50 m (Finnigan, 2000;Raupach et al., 1996), requiring a horizontal filter width of around 5 m. Smaller filters may be required in the lee of hills, around which L s can take smaller values (Finnigan, 2000;Ross, 2008). If gaps are to be simulated explicitly (i.e., openings with O d ≥ L s ), then a horizontal filter width of around O d /8 or smaller is required to adequately resolve the turbulent structures in the gaps. Slightly coarser resolutions than those suggested may suffice for very unstable conditions, but finer filters may be required in very stable cases, because L s decreases with increasing stability. In very stable and very unstable conditions around forests, mixing-layer-type eddies cease to dominate the turbulence structure. In very unstable conditions, the filter width should be chosen to scale with the dominant convective plumes, rather than L s . It is not clear whether the distributed-drag parametrization is appropriate to model forest-atmosphere exchange in very stable conditions, during which the turbulent exchange can be dominated by intermittent bursts rather than mechanically generated eddies (Section 2.6; Mahrt, 2014).
Although LES resolution has increased more slowly than the general semiconductor capacity predicted by Moore's law (Bou-Zeid, 2014), computing capacity is expected to increase with time for the foreseeable future. This raises the question of whether the extra capacity should be applied toward increasing the resolution or the size of the simulated domain. In many fluid applications, the answer to this question is "both." However, the distributed-drag parametrization assumes the averaging volume-the filter width for LES-to be much larger than the individual tree and plant elements so that their presence is accounted for statistically. For mature forests, where the tree trunks can have diameters of a meter or more, decreasing the horizontal filter widths below a few meters makes little physical sense without a formal re-think of the averaging operations. In many cases, the increased computational capacity could be used more productively (a) to increase the simulated ABL height, thereby increasing the total amount of resolved TKE (Grylls et al., 2020), (b) to increase the size of the simulated domain, (c) to include a wider range of physical and chemical processes, or (iv) to account for plant movement. However, there may be instances for which a very fine filter is required, such as to investigate turbulent transfer in small gaps, in the lee of steep hills, or in stable conditions.
Outside of these specific applications, the workhorse spatial resolution of Δ x = Δ y ≈ 5 m, and Δ z ≈ 1 m, means a forest volume of up to tens of cubic meters is accounted for by a single grid cell, with the smallest resolved eddies ≈5 m in diameter. This implies a need for robust SGS schemes if LES simulations are to be useful. Although LES has become popular over the last few decades, the development of SGS schemes has been slow (Moser et al., 2021). The most widely used SGS parametrization is the Smagorinsky scheme, in which closure is achieved using empirical arguments and theory (Smagorinsky, 1963). This scheme assumes the turbulence is isotropic, which is not the case in forests (and in many other ASL applications). Recently, more suitable schemes have been developed by determining the minimum energy dissipation needed to balance the turbu-lence production at SGS scales (Gadde et al., 2021;Rozema et al., 2015), but these schemes have not yet been widely adopted. Improving SGS schemes for forest applications is an important line of future research. Eddies larger than the plant elements lose TKE to heat and fine-scale TKE, which short cuts the usual eddy cascade that is assumed by Kolmogorov's −5/3 law (Finnigan, 2000). Shaw and Patton (2003) show this effect can be accounted for in LES using a TKE cascade term, without the need for extra modifications to account for wakescale kinetic energy transfer. However, the effect of this short cut to fine scales has not been tested in patchy vegetation canopies, with moving plants, or with a view to determining its impact on scalar exchange across the boundary layers of leaves.

Future Challenges
It is an exciting time to investigate forest-atmosphere exchange. There is a growing body of high-quality observations from micrometeorological campaigns (e.g., Butterworth et al., 2021;Patton et al., 2011), long-term ecosystem monitoring sites (e.g., Hicks & Baldocchi, 2020;Oliphant, 2012), and ecosystem-scale FACE experiments (e.g., Drake et al., 2016;Hart et al., 2020;Norby et al., 2006). These observations allow researchers to investigate an increasingly diverse range of physical and biological phenomena in an increasingly diverse range of ecosystems. Numerical models, if deployed correctly, offer a fourth pillar in ecosystem-scale investigations, alongside measurements, theory, and physical models. However, fulfilling the potential of observational networks and new modeling techniques requires a concerted effort across scientific disciplines. Moving beyond idealized cases and moderate atmospheric conditions is not easy. First, the idealizations and simplifications are usually for good reason, for example, to reduce the dimensionality of investigations, or to attempt to reveal canonical processes rather than the idiosyncrasies of individual sites. Second, because many of the experimental and modeling techniques were developed from studies of idealized forests in moderate conditions, attempting to interpret results outside of these scenarios is not always straightforward because the techniques are not at their best (Baldocchi, 2008;Belcher et al., 2012;Brunet, 2020;Stoy et al., 2013;K. Wilson et al., 2002). Third, research into ecosystem-scale forest-atmosphere interactions is motivated by scientific questions across disciplines including forest ecology, fluid mechanics, plant physiology, meteorology and climatology, air and water quality, as well as policy and commercial need. The literature is therefore vast and widely dispersed. In a sense this dispersion is exciting because it suggests that progress may lie at the intersections of scientific disciplines. However, it is also a barrier to progress because results familiar to one scientific community may be unknown in another, even when they relate to important coupled phenomena. We conclude by framing progress and outstanding challenges under three overlapping headings.

Customized Numerical Models of Real Sites
Advances in non-destructive scanning techniques, computing power, and theory are poised to allow researchers to investigate forest-atmosphere exchange at real sites, using models capable of resolving turbulence. These models could be calibrated to study a variety of ecosystem-scale processes, including the wind dispersal of spores, turbulent fluxes of important scalar quantities, or the forest's microclimate. Regarding the structural element of this challenge, TLS allows researchers to include detailed, site-specific morphological detail in their models (Boudreault et al., 2015;Raumonen et al., 2015;D. Wang et al., 2020), and virtual canopy generators can be used to generate realistic forest canopies from a small number of structural variables (Bohrer et al., 2007). Researchers should include a careful description of the forest morphology in their numerical models, particularly when assessing simulated results against observations, or in forests with a high proportion of edge region. Ecosystem-scale exchange models should account for the streamlining of plants in high wind conditions. This can be accounted for efficiently in LES simulations by incorporating the Vogel number B into the parametrization of the aerodynamic drag f i . Plant movement at an ecosystem scale can be incorporated through poroelastic and mechanical parametrizations (Dupont et al., 2010;Pivato et al., 2014) or by the inclusion of a stochastic forcing into the momentum equations at a frequency corresponding to plant movement.
However, two important knowledge gaps need to be narrowed before these customized models can be used to investigate the exchange of scalar quantities at real sites. The first is that there are few public datasets against which nascent models can be tested. The ability of turbulence-resolving models to simulate scalar transport in real landscapes therefore remains essentially unknown. The CHATS observations are the apogee of surface meas-urements in a tree dominated landscape (Dupont & Patton, 2012;Patton et al., 2011) and a valuable resource for model developers (Ma & Liu, 2019). However, the CHATS measurements were made in a homogeneous orchard on level ground and are therefore unsuitable for testing models' ability to handle patchiness, species diversity, and undulating topography. More measurements are needed, especially of scalar quantities such as CO 2 , water vapor, pollutants, and spores in patchy forests. For practical reasons, these measurements may not be on the scale of CHATS. However, existing eddy-covariance and FACE facilities have generated extensive timeseries that may be exploited, especially if combined with experiments of opportunity and targeted observational campaigns around edges, gaps, and hills. For example, further experimental testing will determine the extent to which scalar fluxes reach an approximate equilibrium in patchy landscapes, and whether the chimney effect in the lee of hills is observed in nature as well as in numerical models.
The second gap concerns the soil-forest-atmosphere coupling of scalar quantities. Researchers are beginning to incorporate multi-layer canopy exchange schemes in ecosystem-scale LES models (Ma & Liu, 2019;Patton et al., 2016) and in models designed for slightly larger space and time scales (Bonan et al., 2021;L. Xu et al., 2017;Yan et al., 2020). These schemes include parametrizations of phenomena such as evapotranspiration, soil fluxes, and the photosynthetic uptake of CO 2 that are designed to hold approximately across different ecosystems. A potentially fruitful area of research would be to modify these schemes to reflect observations at individual sites, for example, by adjusting the underlying parametrizations, or through some combination of Bayesian data assimilation and machine learning. A point of caution is that observations in forests tend to be sparse in space (e.g., wind measurements) and sometimes time (e.g., spore releases), and are often measured using proxies or with indirect techniques (e.g., photosynthetic capacity; Gardner et al., 2021). Machine learning techniques should therefore be carefully adapted to account for uncertainty in our knowledge of the underlying physical processes (Geer, 2021).

Connection to Larger Scales
Developments in theory and computing capacity are allowing researchers to begin incorporating ecosystem-scale phenomena into numerical weather prediction and earth system models (Arthur et al., 2019;Bonan et al., 2018Bonan et al., , 2021Shao et al., 2013;Yan et al., 2020). To be computationally feasible, these schemes rely on simplified versions of the forest canopy (Yan et al., 2020), a priori turbulence parameters (Y. Chen et al., 2016), or modifications of MOST with mixing-layer length scales (Bonan et al., 2018). These approximations perform best when the atmosphere is neutral and the surface homogeneous. Further work is needed to determine whether the length scales and approximations can be modified to account for patchy landscapes and a larger range of atmospheric conditions. Robust parametrizations for scalar quantities are especially difficult to find. Missing theoretical links may become apparent through further testing of customized LES models against high-quality field measurements, as outlined above. Another relatively unexplored approach is to reject the assumption that the scalar statistics should always relate to those of the velocity field. For example, understanding the movement of water vapor around at fragmented forest at dusk is a problem that may not yield to a deterministic, drag-based treatment or approximate vertical turbulence profiles. The geometries of velocity and scalar timeseries are often much less scale-dependent than their accompanying physics (Belušić & Mahrt, 2012;D. Chen et al., 2019;Kang, 2015). Turbulent events in these timeseries can also be clustered by their statistical properties, with no assumption of the underlying physical structures Sun et al., 2015). The timeseries of well-chosen turbulence measurements may reveal scale-independent behaviors to parametrize forest-atmosphere exchange in terrain and conditions that are beyond the reach of current approaches. Spectral analysis of long-term eddy-covariance measurements over mixed hardwood forests showed smaller eddies are relatively more efficient in transporting sensible heat than momentum, while larger eddies are relatively more important in transporting momentum (Su et al., 2004). It would be fascinating if this type of analysis were conducted over even longer measurement periods, with a larger diversity of ecosystems.

Challenging Weather and Atmospheric Conditions
Numerical models permit low dimensionality experiments that are difficult to perform in the open atmosphere. A common simplification is that of an isolated forest in neutral atmospheric conditions. This removes all aspects of the atmosphere, other than a velocity field that is recycled using periodic boundary conditions. In nature, however, forests are subject to all sorts of weather, climatic, atmospheric conditions. Atmospheric stability can change quickly in time, such as around dusk, and in space, such as when the forest air space decouples from the surrounding ABL. Simulations of the entire ABL will help explore these effects further (Ma & Liu, 2019;Patton et al., 2016;Yan et al., 2020). Local effects around gaps and areas of uneven heating should not be discounted and, with time, may be possible to resolve in site specific models. However, further research is needed into the effect of submeso motions on forest-atmosphere exchange. These motions are difficult to model because they often manifest as intermittent turbulent bursts, but do not result from mixing-layer-type eddies, and therefore must be generated using mesoscale coupling or synthetic turbulence. Submeso motions are common during the atmospheric conditions that tend to be most problematic for eddy-covariance measurements, such as stable nights, and discounting their effect from models introduces an unwelcome bias. Rainfall is a significant source of momentum into forests, affecting myriad ecological processes as well as the plants' mechanical response to the flow. Raindrop-induced vortices can carry small particles away from plant surfaces and into the turbulent flow. Wet conditions provide a major unexplored area for numerical investigations of forest-atmosphere exchange, especially given that much of the world's forested area is situated in climates where precipitation is common.

A1: Drag Parametrization and Simulated Cases
We simulate flow across small forests, with the presence of the aerial parts of the forests represented using the drag parametrization in Equation 7, = − ( )|̃| <̃> (the overtilde denotes the resolved variables, as in Section 2). We specified a height-averaged equational drag coefficient ̄= 0.2, a value commonly used in previous studies (Table A1). Figure A1 shows the mean vertical profile of the LAD, a(z), which was derived from terrestrial lidar surveys of the BIFoR FACE facility (Hart et al., 2020), using the method in Raumonen et al. (2015), giving PAI ≈ 5. We specified two cases: 1. Case 1-We apply the distributed-drag parametrization over a patch extending 500 m in the streamwise direction, and across the entire domain in the spanwise direction, to simulate flow across a small, isolated forest. 2. Case 2-As for Case 1 but, for each forested cell, ̄( ) taking a random number (uniformly distributed) within the range ̄( ) (Case 1) ±2%. This introduces a small spatial stochastic forcing into the momentum equations solved by the LES model via the drag term in Equation 7. The stochastic forcing does not directly vary with time using this formulation.

A2: Transport Equations
We solved the transport equations using the LES mode of WRF version 3.6.1 (Skamarock et al., 2008). The WRF model solves discretized forms of the spatially averaged momentum equations using the Runge-Kutta time-integration scheme (Wicker & Skamarock, 2002), where ρ is the air density; is the resolved pressure; τ ij is the kinematic mean stress tensor, which represents the SGS stresses; B i is the buoyancy force: = − 3 ′ ∕̄ where ̄ is the potential temperature for hydrostatic balance, and θ′ is the temperature variations with respect to ̄ ; f c is the Coriolis parameter; 3 is the alternating unit tensor; and U g,j is the geostrophic velocity. Equation A1b is closed by parametrizing the SGS stress tensor τ ij as where e SGS is the SGS TKE and c k = 0.10 is a modeling constant. The prognostic equation for the evolution of the term e SGS is, where P r represents the shear-and buoyancy-production terms (Skamarock et al., 2008), C ϵ is the dissipation coefficient (Moeng et al., 2007). The term F is a cascade term, which accounts for additional dissipation of kinetic energy from air-forest interactions at scales smaller than the spatial filter (Shaw & Patton, 2003;Shaw & Schumann, 1992), with = −2̄( )|̃| SGS. (A6)

A3: Drag Parametrization and Simulated Cases
The simulated domain is 1,435 × 1,435 × 1,000 m in the streamwise, spanwise, and vertical directions, respectively, comprising 287 × 287 × 79 grid cells. The horizontal grid resolution is 5 m in each direction, and the vertical resolution is increased from around 0.67 m in the lower half of the forest to around 60 m at the top of the simulated domain. The mean height of the forest h c = 25 m. The upstream edge of the forest is situated approximately 600 m from the inflow edge of the domain. We simulate flow under neutral conditions, with the initial profile of potential temperature θ specified as a constant at 283.15 K at the bottom of the domain (up to z ≈ 475 m) followed by a linear increase to 291.7 K at the top of the domain. We include a dampening layer of z ≈ 300 m at the top of the domain to minimize wave reflection (Nottrott et al., 2014). The geostrophic velocity components above the boundary layer top are set to U g = −6 m s −1 and V g = −9.3 m s −1 , and this specification yields a mean wind speed of 1.6 m s −1 from a flow direction of 343° (approximately a northerly wind) at the crown top (z = h c ). We use the meteorological convention where the x-direction is aligned west-east and the y-direction south-north.
Spin-up was 5 hr, with cyclic boundary conditions for all dynamical variables in both horizontal directions. After the spin-up, we ran the simulations for a further 120 min, taking samples at intervals of 3 s. We time-averaged over the latter 100 min (denoted by T) of these samples (i.e., t 0 = 20 min to t 0 + T = 120 min). This process generates a three-dimensional time series of 2,000 resolved samples in the form ̃( ) . We derived the resolved turbulent statistics using (a) time averages over the sampling period and (b) spatial averages along the y-direction (L y = 1,435 m), over which the turbulent statistics are homogeneous. For each resolved variable, ̃ , the averaging process generates a two-dimensional function, from which we calculate the statistics presented in Figure 10. The resolved fluctuating component of ̃ around Φ is defined as ′ ( ) =̃( ) −Φ( ) with the skewness = ′3 ∕ 3 and kurtosis = ′4 ∕ 4 .

Conflict of Interest
The authors declare no conflicts of interest relevant to this study.

Data Availability Statement
The WRF-LES data used in this paper, and their supporting R scripts, are available at https://github.com/foresteddy/forest_exchange. Figure