Vapor pressure deficit helps explain biogenic volatile organic compound fluxes from the forest floor and canopy of a temperate deciduous forest

Biogenic volatile organic compounds (BVOCs) play critical roles in ecological and earth-system processes. Ecosystem BVOC models rarely include soil and litter fluxes and their accuracy is often challenged by BVOC dynamics during periods of rapid ecosystem change like spring leaf out. We measured BVOC concentrations within the air space of a mixed deciduous forest and used a hybrid Lagrangian/Eulerian canopy transport model to estimate BVOC flux from the forest floor, canopy, and whole ecosystem during spring. Canopy flux measurements were dominated by a large methanol source and small isoprene source during the leaf-out period, consistent with past measurements of leaf ontogeny and theory, and indicative of a BVOC flux situation rarely used in emissions model testing. The contribution of the forest floor to whole-ecosystem BVOC flux is conditional on the compound of interest and is often non-trivial. We created linear models of forest floor, canopy, and whole-ecosystem flux for each study compound and used information criteria-based model selection to find the simplest model with the best fit. Most published BVOC flux models do not include vapor pressure deficit (VPD), but it entered the best canopy, forest floor, and whole-ecosystem BVOC flux model more than any other study variable in the present study. Since VPD is predicted to increase in the future, future studies should investigate how it contributes to BVOC flux through biophysical mechanisms like evaporative demand, leaf temperature and stomatal function.


Introduction
Ecosystems are dynamic sources and sinks of biogenic volatile organic compounds (BVOCs), which play important roles in plant defense (Yuan et al. 2009;Junker and Tholl Communicated by Hermann Heilmeier. 1 3 2013), ecological signaling (Schiestl 2010;Holopainen and Blande 2012), oxidation chemistry, and atmospheric particle formation and growth (Fuentes et al. 2000(Fuentes et al. , 2016Kulmala et al. 2013;Faiola et al. 2014). BVOCs are also important in ecosystem and global carbon budgets. Their total global emission is on the order of 0.76 Pg C year −1 (Sindelarova et al. 2014), nearly a quarter of the size of the net land carbon sink (Friedlingstein et al. 2019), and BVOC flux can comprise up to 10% of photosynthetic carbon uptake at the ecosystem scale (Peñuelas and Llusià 2003). These emissions, however, are dynamic and respond to continued shifts in land use, climate, and atmospheric chemistry (Peñuelas and Staudt 2010;Calfapietra et al. 2013;Hantson et al. 2017). It is thus critical to identify primary drivers contributing to ecosystem BVOC fluxes to better understand BVOC dynamics in a changing world.
Most ecosystem BVOC emissions arise from plant foliage ( Fig. 1) (Guenther 1997), particularly via plant stomata (Fall and Monson 1992), with important additional sources from leaf litter (Leff and Fierer 2008;Gray et al. 2010Gray et al. , 2014Aaltonen et al. 2011). Dynamic source/sink behavior is also observed within the soil itself (Cleveland and Yavitt 1997;Insam and Seewald 2010;Bachy et al. 2016;Tang et al. 2019;Rinnan and Albers 2020). For example, in forests, soil BVOC flux results from belowground dynamics including the functioning of roots (Kreuzwieser and Rennenberg 2013;Gray et al. 2014) and their mycorrhizal associations (Trowbridge et al. 2020), which can subsequently alter the composition of the microbial communities that give rise to BVOCs (McBride et al. 2020). It is unclear if soils and litter make a sufficient contribution to total ecosystem BVOC flux to warrant inclusion in ecosystem models (Asensio et al. 2007), though recent studies indicate that they likely do (e.g. Bachy et al. 2016;Mäki et al. 2019). Aboveground BVOC sinks include reactions with oxidative species within the plant canopy (Carter 1994;Fuentes et al. 2007), dry deposition (Guenther 2015) and, critically, uptake by plant stomata such that stomatal BVOC exchange should be treated as bi-directional (Niinemets et al. 2014) (Fig. 1). Despite our recognition of complex interactions among BVOC emissions drivers, ecosystem BVOC flux models tend to focus on emissions from just the plant canopy, rather than including soils and leaf litter (Baldocchi et al. 1999;Niinemets et al. 2013;Guenther 2015), and they are seldom challenged with observations from a broad set of ecosystem components. This study was aimed at enhancing our understanding of ecosystem BVOC emission modeling with regard to soil and litter processes and model inputs.
Models of ecosystem BVOC exchange are also challenged by plant phenology and seasonality (Holzinger et al. 2006). These include emission bursts during periods of new foliage growth (Aalto et al. 2014 that may be due to Fig. 1 Biological and physical sources, sinks, and transport of biogenic volatile organic compounds (BVOCs) in an idealized plant canopy with leaf area density (LAD, m 2 m −3 ) measurements following Oliphant et al. (2006) and canopy profile measurement heights from the study ecosystem: a mixed deciduous forest in the Morgan Monroe State Forest, Indiana, USA the direct exposure of plant resin to the atmosphere (Eller et al. 2013), other aspects of plant-water relations like xylem refilling (Vanhatalo et al. 2015), or the production or methanol and acetone during leaf ontogeny (MacDonald and Fall 1993a;Nemecek-Marshall et al. 1995). These plant-associated fluxes are controlled by a complex interplay among environmental and biological variables, including intrinsic cellular processes as well as extrinsic factors such as air temperature, solar radiation, and soil moisture (Monson et al. 1995;Guenther et al. 2012). Our knowledge about the interactions among plant processes and the environment, however, is continually increasing and has progressed since the creation of the current generation of ecosystem BVOC emission models.
Model development typically lags behind empirical discovery, making it likely that there are additional intrinsic and extrinsic variables that need to be added to existing models. These variables might include the vapor pressure deficit (VPD) which represents atmospheric demand for water and is critical for plant canopy conductance (Novick et al. 2016), wind speed and canopy roughness, which influence canopy aerodynamic conductance and canopy air-space venting (Bohrer et al. 2009), the diffuse fraction of solar radiation, which can penetrate the plant canopy more efficiently than the direct beam to impact the radiative environment of the subcanopy (Oliphant and Stoy 2018;Moon et al. 2020), and the surface (skin) temperature, which may provide a more accurate description of the temperature at which biological processes occur compared to air temperature (Still et al. 2014;Pau et al. 2018). In addition to the emerging importance of soil and litter processes to accurately predict ecosystem BVOC fluxes, there is a need for studies that take a step back from the current form of emissions models and offer a fresh perspective of the relationships among model logic and existing knowledge of canopy and plant ecophysiological processes.
To study the dynamics of BVOC flux during the leaf-out period and investigate micrometeorological and vegetative controls necessary to explain their dynamics, we measured BVOC concentrations within a mixed deciduous canopy and estimated source and sink areas and fluxes using a hybrid Lagrangian/Eulerian canopy transport model. The transport model incorporated BVOC concentration measurements from a unique logarithmic canopy profiling system with a higher density of observations from lower in the canopy airspace, compared to past studies, which allowed us to isolate the contribution of the forest floor to ecosystem-scale BVOC exchange. The biotic diversity of species and traits in mixed forests creates additional challenges to understanding controls over sources and sinks of BVOCs (Kaharabata et al. 1999) and we included analyses of wind direction to understand if flux source area was important to describe efflux in the study forest. Finally, models of BVOC flux are increasingly cognizant of compounds that make a minor contribution to mass flux (Guenther et al. 2012) but may play disproportionate roles in ecological interactions and atmospheric chemistry (Goldstein and Galbally 2007;Clavijo McCormick et al. 2014). We took care to not exclude an analysis of these minor compounds by exploring their relationship with BVOC compounds that comprise a larger proportion of total BVOC flux.

Study site
Measurements were made in the Morgan-Monroe State Forest in south-central Indiana, USA, at 39°190′ N, 86°250′ W on and around a long-running Ameriflux eddy covariance tower (site code US-MMS, Schmid 2000). The study site is an approximately 80-year-old mixed hardwood forest with trees that associate with ectomycorrhizal fungi (ECM) including shagbark and pignut hickory (Carya ovata and C. glabra), red and white oak (Quercus rubra and Q. alba), and American beech (Fagus grandifolia), and those that associate with arbuscular mycorrhizae (AM) including tulip poplar (Liriodendron tulipifera), sugar maple (Acer saccharum), and sassafras (Sassafras albidum). Seedlings and the isoprene-emitting spicebush (Lindera benzoin) are found in the understory. The mean canopy height in the vicinity of the tower was 27 m when measurements were made. Please see Brzostek et al. (2015) for additional site details.

Micrometeorological measurements
The US-MMS eddy covariance tower includes a full suite of micrometeorological measurements including eddy covariance systems comprising CSAT3 (Campbell Scientific, Logan, UT) sonic anemometers and LI-7000 closed path infrared gas analyzers (LICOR, Inc., Lincoln, NE) at 46 m, 34 m, and 2 m above the forest floor; HMP35C air temperature/relative humidity measurements (Vaisala, Vantaa, Finland) at the same heights; a CNR-1 four-component radiometer (Kipp and Zonen, Delft, The Netherlands) at 46 m; and CS615/6 (Campbell Scientific) soil moisture sensors and soil temperature measurements at 10 cm depths to bedrock. We use direct and diffuse photosynthetically active photon flux density (PPFD), measured by a BF3 sunshine sensor (Delta-T Devices, Cambridge, UK) to describe BVOC responses to the light environment.

Proton transfer reaction-mass spectroscopy
BVOC measurements were made using a proton transfer reaction-mass spectrometer (PTR-MS; Ionicon Analytic, Innsbruck, Austria) housed in an air-conditioned research building adjacent to the US-MMS tower. We selectively analyzed 49 different mass/charge (m/z) ratios following previous work on compounds identified in ambient air (Gouw and Warneke 2007;Blake et al. 2009;Ellis and Mayhew 2013) (Table S1). We focus on key compounds for which we have calibration standards, namely formaldehyde, methanol, acetonitrile, acetaldehyde, acetone, dimethyl sulfide (DMS), isoprene, methyl vinyl ketone (MVK), methacrolein, methyl ethyl ketone (MEK), benzene, toluene, C8 aromatics, C9 aromatics, and monoterpenes (Table S1). Concentrations of the fragments of isoprene and monoterpene (m/z 41 and m/z 81) were calculated using a theoretical transmission curve created for the PTR-MS instrument. The concentrations of the fragments were then added to the concentration measurements of the compounds calculated using calibration. Calibrations were performed prior to each measurement campaign using a multi-component calibration mix (± 5% confirmed by the manufacturer using GC-MS) stored in nitrogen gas (Apel-Riemer Environmental, Inc.). We studied the time series of the normalized mass-to-charge (m/z) ratios of compounds for which calibration standards were not available to see if they follow similar patterns to compounds with calibration standards.
The PTR-MS drift tube pressure (p) was set at 2.1 mbar and temperature (T) to 333.15 K with a drift field of 600 V. The parent ion signal was set to ~ 1 × 10 6 counts per second. The protonated oxygen to primary ion count (O 2 + :H 3 O + ) ratio was < 3.5%. The volume mixing ratio (VMR, ppbv) of calibrated compounds was calculated by dividing its m/z by its calibration coefficient calculated using a linear fit. To calculate concentrations using the transmission factor approach, the volume mixing ratio VMR of each compound was calculated following (Ellis and Mayhew 2013): where i MH + is the protonated compound of interest divided by its transmission factor, i H 3 O + is the primary ion count divided by its transmission factor and multiplied by 500, k is the rate coefficient, t is reaction time, and the total number density of the gas in the drift tube, N d , is where p is in Pa, A N is Avogadro's number, and R is the ideal gas law constant (8.3145 m 3 Pa mol −1 K −1 ). We note that this approach incurs uncertainty on the order of 2% using 500 rather than 489.56 for the ion count multiplier. (1)

Canopy profile measurements
The canopy BVOC profile measurement campaign took place from May 8 to May 25, 2015 with a 77-min missing measurement period during maintenance on May 12. During sampling, the inlet of the PTR-MS was attached to the canopy profiling system with measurements at 0.25 m, 0.5 m, 1 m, 2 m, 4 m, 8 m, 16 m, and 32 m (i.e. 2 −2 to 2 5 m). Measurements were made every 10.45 s. The first three and last PTR-MS measurements at each level were discarded to ensure that the sample air originated from the measurement height at the time of measurement. The remaining observations at each level were then used to calculate a mean concentration of each of the study compounds.

Canopy flux estimation
Estimating the net flux of a scalar (here, BVOCs) from canopy profile observations requires a model that can compute spatially distributed scalar sources from time-distributed concentration measurements, known as the 'inverse problem' (Raupach 1989a). The challenge is to infer the source and sink (S) dynamics of a scalar as a function of height (z) from measurements of the scalar concentration (c) at different heights in the canopy across time (t) which are coupled by the principle of continuity and the scalar concentration budget equation (see Monson and Baldocchi 2014). The time-and horizontally averaged steady-state scalar conservation equation for planar homogeneous high Reynolds and Peclet numbers flow (i.e. neglecting molecular diffusion, (Finnigan 1985;Raupach 1988) is given by where w is vertical velocity, primes are fluctuations from time averages (represented by overbars), and the angle bracket denotes horizontal averaging as discussed in Raupach and Shaw (1982), such that ⟨ w ′ c ′ ⟩ represents the turbulent vertical flux of a scalar. Early approaches used 'K-theory', which assumes that scalar fluxes and, consequently, S (Eq. 3), are related to scalar gradients through an eddy diffusivity coefficient within the plant canopy, which can vary with height. However, the dynamics of turbulence within plant canopies, including counter-gradient fluxes, often result in a situation where K-theory becomes unreliable (Denmead and Bradley 1985). Whereas some approaches have improved models for eddy diffusivity to account for these shortcomings (Freire et al. 2017), here we adopt the 'hybrid' approach of Siqueira et al. (2000) that combines the two dominant approaches designed to overcome the limitations of K-theory: inverse Lagrangian localized near-field (LNF) theory (Raupach 1989a, b) and high-order Eulerian closure models (Katul and Albertson 1999). In brief, the hybrid approach acknowledges that LNF theory assumes normally distributed vertical velocity statistics, but also recognizes that non-Gaussian ejection-sweep cycles frequently drive mass transport within tall plant canopies ( Fig. 1), and that pure Eulerian formulations are sensitive to uncertainties in the measurement of c due to a limited scope for replication in the inversion of S from c (Siqueira et al. 2000). To overcome these limitations, Siqueira et al. (2000) proposed a model that retains the concept of the inverse Lagrangian approach, which grants the model the desired robustness to measurement errors but using the Eulerian frame of reference to estimate the elements of the dispersion matrix, thus mechanistically accounting for vertical-velocity skewness effects in scalar transport.
The procedure consists of dividing the canopy into multiple layers (Fig. 1). The task then becomes to find the combination of the canopy layer source and sink strengths (i.e. S) that best recover the measured mean scalar concentration profile. In an inverse sense, the dispersion matrix is constructed by dividing the canopy into layers of unit source strength and then calculating the expected mean concentration profiles from each layer individually using the scalarflux budget equation (Eq. 3), which is re-arranged to provide a differential equation for scalar concentration (see Supplemental Material). Neglecting buoyancy, scalar drag and waving source terms, the time and horizontally averaged steady-state equation can be written as a second-order ordinary differential equation (ODE) for scalar concentration as a function of height (Siqueira et al. 2000): where c is the scalar concentration function of height. The coefficients M 1,2,3 are functions of scalar flux and turbulent velocity statistics. The profiles of velocity statistics, when normalized by friction velocity (u*) at a reference height, become a function only of leaf-area-density vertical distribution (LAD), provided for MMSF by Oliphant et al. (2006) (Fig. 1) and can be obtained by a second-order closure model for turbulent flow. Here, we used the model described in Siqueira et al. (2012), who explicitly solved an equation for turbulent-kinetic-energy dissipation, required for the M-coefficients.
Equation (4) is not closed because it requires the scalar fluxes, which can be obtained by Eq. (3) given source distribution is known. With a prescribed source, as is the case for dispersion-matrix construction, the scalar budget equation (Eq. 3) can be integrated to provide the flux used in the M-coefficients. Furthermore, with proper boundary condition, Eq. (4) can be numerically solved (here we used a finite-volume technique) to give the vertical profile of < − c>. Next, the dispersion matrix is computed from where < − c> i represents the concentration at the measurement height i (i = 1,2,…,n) resulting from the source layer j (j = 1,2,…,m), calculated using Eqs. (3) and (4), and − C R is the concentration at a reference height. We adopted the last measurement height (32 m) as the reference, which was used as the boundary condition for Eq. (4). D ij are the elements of the (n by m) dispersion matrix, s is an assumed unitary source strength, and Δz j is the source layer thickness. Once the dispersion matrix is determined, the source strengths S j can be readily computed if m = n: where − C i,m , contrary to (5), are the time-averaged concentrations at height i. However, this would make estimated source sensitive to measurement errors (Raupach 1989a). To avoid such instability, redundant concentration measurements are necessary (i.e. n > m), such that the system becomes overdetermined. As shown by Raupach (1989a), such redundancy reduces (6) to a regression problem with the source strengths calculated by a least-squares approach: where and The above regression procedure addresses the limitedsampling problem inherent in pure Eulerian frameworks. However, such an approach, when applied to limited sample size, retains high variance among redundant source estimates. To improve the assessment of retrieved source values and their associated variances, an additional smoothness constraint was imposed on (7) and (8) using the Weighted Measures of Length procedure (Menke 2018;Siqueira et al. 2000).

Modeling analysis
We created linear models of the forest floor, canopy, and ecosystem flux of each studied BVOC compound (Table 1) as a function of multiple micrometeorological variables and the surface-atmosphere flux of carbon, water, and heat, fit using maximum likelihood. We discriminated amongst the models by calculating the Akaike's Information Criterion for each with the assistance of the dredge function in the MuMIn package (Bartoń 2020) using R version 4.0.0 (R Core Team 2020) and selecting the model with the lowest AIC as the most parsimonious. The goal of the modeling analysis is to determine which variables should not be excluded for understanding micrometeorological and canopy controls over BVOC flux rather than creating models of the flux of 14 compounds from three sources during the leaf-out period of a mixed deciduous forest. The micrometeorological variables used to create the modelsin addition to those listed above-include below canopy and diffuse PPFD (as a surrogate for incident solar radiation), soil temperature and moisture, VPD, wind speed and direction as a surrogate for source area within the diverse forest canopy, and canopy and forest floor surface (skin) temperature calculated from outgoing longwave radiation measurements using the Stefan-Boltzmann Law and assuming a canopy and forest floor emissivity of 0.98 (Jin and Liang 2006).

BVOC concentrations within the canopy
Methanol had the highest mean (± SD) concentration across all measurement heights at 12.7 ± 6.5 ppbv (Table 1). Within-canopy mean BVOC concentrations of the key compounds methanol, isoprene, and monoterpenes were relatively high during the earlier part of the measurement period (Fig. 4), then declined and increased again on May 15 during a warm period after a small rain event. Canopy BVOC concentrations rose again during the warm and sunny period of May 23-25. Notably, methanol, isoprene, and monoterpene concentrations were often higher in the overstory at 16 m during May 8-9 and May 13 and isoprene and monoterpene concentrations were often higher near the soil surface and subcanopy at 0.25 m and 0.5 m from May 15-25 (Fig. 4), implying different sinks and sources within the canopy volume.
The time series for most key compounds exhibited characteristic concentration profiles as a function of time of day and canopy height with elevated concentrations in the afternoon within the canopy and near the forest floor (Fig. 5) except methanol, which had an early morning peak throughout the canopy volume, on average, and formaldehyde which had a peak in the lower forest canopy in the afternoon. The mass-to-charge ratios of most of the 49 compounds studied also exhibited a characteristic pattern with higher values within and above the canopy and in the afternoon ( Figure S2) noting that many compounds were close to the signal-to-noise ratios as evidenced by Table 1 The mean and standard deviation of studied biogenic volatile organic compound (BVOC) concentrations across all measurement heights ( Fig. 1) Figure S2).

Modeled BVOC fluxes, sources, and sinks
Models of the second and third moments of the vertical velocity normalized by u* matched observations well on average ( Figure S3). Modeled BVOC concentrations likewise tended to match measurements well with the exception of the lower canopy where modeled values were frequently lower than measurements (Fig. 6) due likely in part to challenges in modeling turbulence near surfaces (the boundary layer effect) and challenges in redistributing parcels toward the surface in Lagrangian approaches. Forest floor BVOC source strength may be underestimated and/or sink strength overestimated as a consequence and results are subject to this uncertainty. Most BVOC compounds exhibited a relatively large net efflux to the atmosphere from the plant canopy from the beginning of the measurement period starting at mid-day on May 8 until mid-day on May 9 ( Fig. 7) that resulted from (relatively) high BVOC concentrations in the overstory at 16 m during the first day of measurements (e.g. Fig. 4). There was a subsequent net uptake of most studied BVOC . a Air temperature measured at 46 m (T a ), soil temperature measured at 10 cm depth (T s ) and above-canopy precipitation (P). b Photosynthetically active photon flux density (PPFD). c Vapor pressure deficit (VPD). d Wind direction (WD) measured above and within the plant canopy compounds from May 9 until May 12 during a period when mid-day temperature decreased from nearly 28 °C to less than 16 °C with cloudier conditions (Figs. 2, 7). This was followed by another net efflux of most BVOC compounds from mid-day on May 12 until mid-day on May 13 during a sunnier period (Figs. 2, 7). Afterward, net BVOC flux from the overstory tended to be minor for some compounds (e.g. methanol) or exhibit net uptake for others (e.g. isoprene). The forest floor was a net cumulative sink of BVOCs across most of the measurement period, which tended to buffer net canopy fluxes such that the net ecosystem source was smaller than net canopy source alone (Fig. 7), keeping in mind that the forest floor sink strength is likely underestimated (Fig. 6). There was an increase in forest floor BVOC efflux during the last 2 days of measurements from midday on May 23 until mid-day on May 25 (Fig. 7) during a warmer period (Fig. 2a) when subcanopy and above-canopy wind direction was decoupled (Fig. 2d) as the forest canopy approached closure (Fig. 3). As a result of this similar behavior across time, the ecosystem fluxes of all study compounds were significantly related to each other and the flux of one compound explained up to 81% of the variability of other compounds (Figures S4 and S5). The canopy overstory was a net source of BVOCs and the canopy air space below it was on average a net sink with the exception of the understory vegetation which was a notable source of BVOCs, especially toward the latter part of the measurement period as demonstrated for the key compounds methanol, isoprene, and monoterpenes (Fig. 8).

Ecosystem models
Leaf area index (LAI) and VPD entered the most parsimonious ecosystem, canopy, and forest floor BVOC flux model on 39 of 42 and 41 of 42 instances, respectively (Table 2) and had the strongest correlation with the key study species methanol, isoprene, and monoterpenes (Fig. 9). Air (canopy) temperature entered the most parsimonious model on 31 (30) of 42 instances and soil temperature and the diffuse fraction of photosynthetically active radiation entered the most parsimonious model for forest floor flux on 12 of 14 instances. Below-canopy radiation was never an important input, but wind speed and direction often were.

Ecosystem BVOC efflux
We observed a relatively large BVOC efflux from the canopy during the early part of the May measurement period (May 8-9), especially at 16 m (Figs. 4,8). This height roughly aligns with the canopy layer of maximum leaf-area density located vertically near the middle of the foliated portion of the canopy (Fig. 1). The observation of this early-spring efflux occurs during the period of rapid leaf expansion, in this case during the time that leaves expanded from approximately 70% to 95% of full expansion area (Fig. 3). Methanol fluxes were at least an order of magnitude higher than all other BVOC emissions, including isoprene, which normally dominates emissions in eastern US deciduous forests (Geron et al. 1994;Isebrands et al. 1999). It is likely that the high methanol fluxes that we observed, and the fact that isoprene emission rates were much lower, are associated with the physiological maturation of leaves. Methanol is known to be formed at high rates during leaf expansion as a product of pectin demethylation during cell-wall loosening (Levy and Staehelin 1992;MacDonald and Fall 1993b;Galbally and Kirstine 2002). High methanol emissions have been observed during the early spring leaf-out period from other mixed forests in the north-central US (Karl et al. 2002;McKinney et al. 2011) and boreal forests (Aalto et al. 2014;Schallhart et al. 2018). At the same time, isoprene is known to be emitted at low to negligible rates early during leaf expansion and is only fully activated as leaves near full expansion (Grinspoon et al. 1993), or after treatment by relatively high accumulated temperature (Monson et al. 1994). Fluxes and atmospheric concentrations of monoterpenes were relatively low (Table 1; Fig. 7) likely due to the lack of coniferous species in the forest canopy. The combination of relatively high methanol fluxes and relatively low terpenoid fluxes represents a chemical-flux landscape seldom studied with regard to model testing and provided us with an opportunity to challenge the coupling between modeled emissions and climatic drivers in a novel ecosystem context (see also Aalto et al. 2014).
The estimated rate of canopy isoprene emissions at its peak in early May was approximately 10% of the rate previously measured in eastern US forests (Goldstein et al. 1998;Baldocchi et al. 1999). Of course, comparisons of isoprene emissions rates among sites will depend on the fraction of trees at each site that emit isoprene. However, the Morgan-Monroe forest has a relatively high representation of oaks, which are high isoprene emitters, and there is no reason to suspect that such low isoprene emission rates, compared to other sites, are due to forest species composition. Rather, it is more consistent with suppressed emission rates due to leaf ontogenetic effects, as described above. Even with low basal emission capacities, however, isoprene emissions responded to changes in seasonal weather conditions. Isoprene emissions reached relatively high rates during the warm, cloudless period of May 8-12, consistent with its high sensitivity to PPFD and leaf temperature, but the cold, cloudy periods between May 12-16 and again between May 19-22, appear to have caused a persistent decrease, consistent with past studies that have shown a close coupling of isoprene emissions to prevailing weather periods and photosynthesis, and a lag in recovery to high emission rates following short periods  (Sharkey et al. 1999). On May 12, during a period of relatively high solar radiation, the prevailing wind direction switched from its normal SSW flow to become progressively more northerly, causing the profiling system to measure different areas of the forest (Fig. 2b, d). This shift may have contributed to the noticeable dip in emissions of all BVOCs, but especially for isoprene and methanol, on May 12.
BVOC sources dominated the observed total canopy flux, as forest floor sinks were small during this part of the growing season (Fig. 7). Overall, the emission profile for this springtime campaign is skewed toward high methanol emissions, but still responsive to climate variation, especially temperature, with regard to terpenoid emissions. This creates a novel set of data for testing emissions models because there is evidence of clear responses to the conventional climate drivers of temperature and PPFD, but with the added early-season condition of high methanol emissions due to phenological drivers.

Forest floor BVOC efflux
The forest floor was a net sink of most study BVOCs during most of the measurement period (Fig. 7) in agreement with the results of soil BVOC flux measurements from MMSF during the early growing season of the previous year (Trowbridge et al. 2020) and noting again that the forest floor sink strength from profile measurements should be interpreted as an underestimate of the true flux (Fig. 6). These observations align with numerous recent studies demonstrating that soils are often net sinks for BVOCs (Rinnan and Albers 2020;Trowbridge et al. 2020). Importantly, the role of the forest floor as a BVOC sink in the absence of fresh litter inputs during the leaf-out period buffered canopy BVOC flux during most of the measurement period such that wholeecosystem BVOC flux was lower than canopy BVOC flux (Fig. 7), but this effect varied by compound. Monoterpene flux from the forest floor was trivial (a fraction of a mmol m −2 over the study period) such that canopy and ecosystem  (Table S1) as a function of height in the canopy during the May 2015 sampling period at the Morgan Monroe State Forest eddy covariance tower (US-MMS) effluxes were nearly identical (Fig. 7) but the forest floor and canopy flux of compounds like acetonitrile were of similar orders of magnitude such that including forest floor flux is necessary to describe whole-ecosystem flux even if the magnitudes of the flux of these compounds are relatively small (Fig. 7). These results suggest that the inclusion of the forest floor to whole-ecosystem BVOC flux is conditional on the compound of interest during the leaf out period. Notably, many minor compounds tended to be highly correlated to two of the calibrated compounds, DMS (Whelan and Rhew 2016) and acetaldehyde (Karl et al. 2002) (Table S1), that had non-trivial contributions of forest floor BVOC flux to ecosystem BVOC flux across most of the measurement period (Fig. 7) suggesting that forest floor fluxes of minor compounds may likewise be a non-trivial proportion of their whole-ecosystem flux.

Ecosystem modeling
Air temperature and PPFD (as a surrogate of shortwave radiation) frequently entered the most parsimonious model of BVOC fluxes as anticipated (Table 2) and there was little evidence that alternate measurements of temperature or PPFD (e.g., skin temperature or diffuse fraction) represented an improvement (Table 1): air and radiometric canopy temperatures entered the most parsimonious ecosystem, canopy, and forest floor models on the same number of instances but the latter is a more difficult measurement to make. These observations suggest-at least for the study ecosystem and measurement period-that the variability in canopy and forest floor temperature as well as below-canopy and diffuse radiation provide little new information to canopy and ecosystem BVOC models than simple air temperature and incident radiation at or near the top of the canopy (Guenther et al. 2006;Arneth et al. 2011). Likewise, radiometric forest floor temperature did not represent an improvement over soil temperature for forest floor BVOC flux modeling.
Atmospheric VPD, on the other hand, consistently entered the linear model with the lowest AIC score for canopy, forest floor, and whole-ecosystem BVOC flux models, and entered these models more than any other variable, including air temperature and PPFD (Table 2). This is perhaps an unexpected result given that nearly all flux modeling over the past three decades for isoprene and monoterpenes-the dominant compounds emitted from most forested ecosystemshas been founded on temperature and light as the dominant driving variables, and with good physiological justification (Monson et al. 2012). It is possible that the importance of VPD arises because of the diversity of BVOCs that we analyzed and the fact that several of them have low Henry's Law volatility coefficients that renders their flux susceptible to stomatal control. For the case of methanol, the BVOC emitted at highest rates during the spring campaign, it has long been known that leaf emission rates are determined by  (Nemecek-Marshall et al. 1995), or a combination of stomatal conductance and leaf temperature (Harley et al. 2007). The tendency for a BVOC compound to be controlled, or not, by stomatal conductance dynamics was explained using chemical theory by Niinemets and Reichstein (2003). Oxygenated BVOCs, such as methanol, acetone, acetaldehyde, MVK, formaldehyde and acetonitrile, preferentially partition into the liquid phase of the leaf and are emitted in a pattern similar to that of water molecules, with significant modification by stomatal control. Hydrocarbon compounds, such as the terpenoids, partition preferentially into the gas phase of the leaf and are not susceptible to stomatal control during steady-state emissions (also see Fall and Monson 1992). It is also worth noting that the canopy flux of oxygenated compounds, like methanol, is also highly susceptible to uptake into moisture films on canopy surfaces (Laffineur et al. 2012). This can cause the canopy to be a sink for these compounds, especially during the early morning and after rain events, when VPD would also be low. Between stomatal control over emissions at high VPD and uptake to the canopy at low VPD, a correlation between net emission rates and VPD dynamics in the modeling is likely explained (Fig. 9). Thus, the tendency for VPD to control canopy BVOC emissions during the spring might result from it falling into the model for many of the oxygenated compounds that were observed. It is surprising, however, that VPD emerged as a significant control in the modeling of nearly all compounds from all sources (Table 2). This means that it also influenced significant control over dynamics in several of the relatively hydrophobic terpenoids. It is not clear at this time as to how such control occurs. It could be due to the fact that there is a high degree of correlation between VPD and air temperature (Fig. 9), and air temperature exerts such strong control over hydrocarbon emissions (due to their low boiling points). It is possible that the cross-correlation between VPD and temperature is causing VPD to appear as important in the Akaike's Information Criterion analysis. Ecosystem BVOC models tend to use soil moisture to simulate plant water stress (Guenther et al. 2012) and SWC entered many of the most parsimonious BVOC flux models even though it was not lower than 33% during the measurement period and therefore not likely to be limiting (Rodriguez-Iturbe et al. 2001) but it did not enter models as frequently as VPD. Our results suggest that VPD is a logical variable to add to BVOC flux models and we recommend additional experiments to explore its role in BVOC flux at the ecosystem scale.
It is also important to note that wind speed and direction consistently entered models with the lowest AIC values, implying that the source area of the sampled air mass in the diverse study forest is important for describing BVOC flux (Table 2), as anticipated given the importance of source area for BVOC flux measurements (Guenther et al. 1996). LAI was also an important variable modeling BVOC flux as anticipated given its critical role in existing flux models (Guenther et al. 2006). It also entered all models of forest floor BVOC flux suggesting that it may be an effective surrogate of whole-ecosystem BVOC dynamics due to simultaneous belowground autotrophic activity. It is also important to note that PPFD DIF entered the most parsimonious forest floor flux model in most instances, but below-canopy radiation itself did not. These results are consistent with the notion that BVOC flux associated with photodegradationknown to be important to litter decomposition (Austin and Vivanco 2006)-played little role in forest floor BVOC flux but that the light environment below the canopy itself did (Moon et al. 2020). Fluxes of the study BVOC compounds were often highly correlated with each other (Figs. 9 and S4, Table S1), as has been found in multiple other studies (Schade and Goldstein 2001), further lending confidence to the notion that BVOCs can be modeled categorically (Guenther et al. 2012). As a whole, our modeling results point to the importance of the canopy light environment and evaporative demand for controlling BVOC flux in addition to the key variables included in models.

Conclusions
We coupled BVOC flux estimates from canopy profile observations and a canopy transport model. Modeled BVOC concentration tended to fit observations well with the exception of the lower canopy layers, suggesting that the modeled forest floor BVOC sink may be underestimated. Observations demonstrate that the addition of VPD may be a logical approach for further improving BVOC model fit and that the contribution of the forest floor to whole ecosystem BVOC flux was either trivial or non-trivial depending on the compound of interest. The fluxes and concentration time series of many compounds were highly correlated, further lending strength to the idea that they can be modeled categorically. Future research should further explore the mechanisms by which VPD controls ecosystem BVOC flux and how the forest floor and canopy combine to create whole-ecosystem BVOC fluxes.

Table 2
The number of times that biological and micrometeorological variables entered the linear model with the lowest AIC for models of the flux of each study biogenic volatile organic compound (BVOC ,  Table 1) FC eddy covariance-measured carbon dioxide flux, H sensible heat flux, LAI leaf area index, LE latent heat flux, PPFD photosynthetically active photon flux density, PPFD BC below-canopy PPFD, PPFD DIF diffuse PPFD, SWC soil water content, T A air temperature, T CAN radiometric canopy temperature, T FF radiometric forest floor temperature, T S soil temperature, VPD vapor pressure deficit, WD wind direction, WS wind speed  . 9 The correlation between the ecosystem-scale flux of the three main BVOC compounds studied here-methanol, isoprene, and monoterpenesbetween themselves, ecosystemscale carbon, water, and energy fluxes, and micrometeorological variables visualized using 'corrplot' .
Colors and ellipsoid shapes correspond to correlation coefficients, shown in the color bar. FC eddy covariance-measured carbon dioxide flux, H sensible heat flux, LE latent heat flux, LAI leaf area index, PPFD photosynthetically active photon flux density, PPFD BC belowcanopy PPFD, PPFD DIF diffuse PPFD, SWC soil water content, T A air temperature, T CAN radiometric canopy temperature, T FF radiometric forest floor temperature, T S soil temperature, VPD vapor pressure deficit, WD wind direction, WS wind speed