A model-based analysis of foliar NOx deposition

Foliar deposition of NO2 removes a large fraction of the global soil-emitted NOx. Understanding the mechanisms of NOx foliar loss is important for constraining surface ozone, NOx mixing ratios, and assessing the impacts of nitrogen inputs to ecosystems. We have constructed a 1D multi-box model with representations of chemistry and vertical transport to evaluate the impact of leaf-level processes on canopy-scale concentrations, lifetimes, and canopy fluxes of NOx. Our model is able to closely replicate 10 canopy fluxes and above-canopy NOx daytime mixing ratios observed during two field campaigns, one in a western Sierra Nevada pine forest (BEARPEX-2009) and the other a northern Michigan mixed hardwood forest (UMBS-2012). We present a conceptual argument for the importance of NO2 dry deposition and demonstrate that NO2 deposition can provide a mechanistic explanation for the canopy reduction of NOx. We show that foliar deposition can explain observations suggesting as much as ~60% of soil-emitted NOx is removed within forest canopies. Stomatal conductances greater than 0.1 cm s result 15 in modelled canopy reduction factors in the range of those used in global models, reconciling inferences of canopy NOx reduction with leaf-level deposition processes. We show that incorporating parameterizations for vapor pressure deficit and soil water potential has a substantial impact on predicted NO2 deposition in our model, with the percent of soil NOx removed within one canopy increasing by ~15% in wet conditions compared to dry conditions. NO2 foliar deposition was also found to have a significant impact on ozone and nitrogen budgets under both high and low NOx conditions. 20


Introduction
The chemistry of nitrogen oxides (NOx ≡ NO + NO2) has a large impact on the oxidative capacity of the atmosphere and the budget of global surface ozone (Crutzen, 1979). NOx is primarily removed from the atmosphere by chemical reactions to form nitric acid, alkyl nitrates, and peroxynitrates, and by dry deposition of NO2 (Crutzen, 1979;Jacob and Wofsy, 1990;Romer et al. 2016). The chemical loss pathways of NOx have been extensively studied, but the physical loss of NO2 to dry deposition 25 remains much more uncertain. Globally, foliar deposition of NO2 removes 20-50% of soil-emitted NO (Jacob and Wofsy,1990;Yienger and Levy, 1995), and constrains near-surface NOx concentrations and input to ecosystems (Hardacre et al. 2015). Understanding the processes that control this removal of NOx by the biosphere is important for predicting anthropogenic surface ozone and understanding flows in the nitrogen cycle.
Reactive nitrogen oxides also serve as an important nutrient in ecosystems. Exchange processes cycle nitrogen between the biosphere and atmosphere, influencing the availability of nitrogen to ecosystems (Townsend et al., 1996;Holland et al., 1997;Galloway et al., 2004;Holland et al., 2005). Deposition of atmospheric reactive nitrogen species can fertilize ecosystems with limited nitrogen availability (Ammann et al., 1995;Townsend et al., 1996;Williams et al., 1996;Holland et al., 1997;Galloway et al., 2004;Teklemariam and Sparks, 2006). Although nitrogen is often the limiting nutrient for plant 5 growth (Oren et al., 2001;Galloway et al., 2004), anthropogenic activities have in some cases caused an excess loading of nitrogen to ecosystems, leading to dehydration, chlorosis, soil acidification, and a decline in productivity (Vitousek et al.,1997;Fenn et al., 1998;Galloway et al., 2004).
The current understanding of the exchange of nitrogen oxides between the atmosphere and biosphere remains incomplete. Despite the importance of dry deposition processes, they are among the most uncertain and poorly constrained 10 aspects of atmosphere-biosphere nitrogen exchange and the tropospheric budgets of O3 and NOx (Wild, 2007;Min et al., 2014;Hardacre et al., 2015). This uncertainty arises from the complex dependence of dry deposition processes on surface cover, meteorology, seasonal changes in leaf area index (LAI), species of vegetation, and the chemical species carrying odd-N.
Developing a mechanistic understanding of dry deposition of NO2 has largely depended on inferences from scarce long-term field observation data and a limited number of laboratory studies on the effects of environmental factors on deposition at the 15 leaf-level. This understanding is represented by a deposition velocity, Vd.
The deposition velocity of NO2 to vegetation is largely regulated by stomatal conductance (Johansson, 1987;Thoene et al., 1991;Rondon and Granat, 1994;Teklemariam and Sparks, 2006;Chaparro-Suarez et al., 2011;Breuninger et al., 2012;Delaria et al., 2018), which varies with tree species, photosynthetically active radiation (PAR), vapor pressure deficit (VPD), temperature (T), soil water potential (SWP) and seasonality of leaf phenology (Emberson et al., 2000;Zhang et al., 2003;20 Altimir et al., 2004;Hardacre et al., 2015;Kavassalis and Murphy, 2017). Many global scale chemical transport models (Jacob and Wofsy, 1990;Wang and Leuning, 1998;Ganzeveld et al., 2002) parameterize Vd using the resistance in-series approach similar to that developed by Wesely (1989). These treatments are heavily parameterized, leading to a large degree of uncertainty, and do not account for the effects of VPD, SWP, CO2 mixing ratio, or other factors known to influence stomatal conductance (Hardacre et al., 2015). NO2 deposition remains even more uncertain than deposition of O3, where stomatal 25 response has been shown to be the primary regulator of foliar deposition and mesophyllic resistance to deposition is negligible.
Observations from leaf-level laboratory studies suggest the deposition of NO2 is also controlled by stomatal aperture (Hanson and Lindberg, 1991;Rondon and Granat, 1994;Hereid and Monson, 2001;Teklemariam and Sparks, 2006;Pape et al., 2008;Chaparro-Suarez et al., 2011;Breuninger et al., 2012;Delaria et al., 2018), however, reactions in the mesophyll may also be important for controlling the deposition velocity of NO2 (Teklemariam and Sparks, 2006;Breuninger et al., 2012). A failure 30 to consider the effects of relevant meteorology on stomatal conductance, as well as our deficient understanding of mesophyllic resistances and the diversity of ecosystem responses, severely limits our ability to understand dry-deposition processes and how they will be affected by feedbacks from changes in climate, land use, and air pollution.
The importance of these considerations has recently been illustrated by Kavassalis and Murphy (2017), who found a significant correlation between VPD and ozone loss, and demonstrated that modeling using VPD-dependent parameterizations of deposition better predicted the correlation they observed. Previous work by Altimir el al. (2004) and Gunderson et al. (2002) have described the effects of VPD and other environmental parameters on the stomatal conductance to O3 of Pinus sylvestris and Liquidambar styraciflua, respectively. More recent models, like the DO3SE model for estimating stomatal conductance 5 to predict ozone deposition velocities, fluxes and damage to plants, incorporate the effects of VPD and SWP on stomatal conductance, but no similar model exists for assessing these effects on NOx deposition. The DO3SE has successfully been implemented in the European Monitoring and Evaluation Program (EMEP) regional model (2012). Modelling studies by Buker et al. (2007) and Emberson et al. (2000) have also demonstrated the success of regional-scale parameterizations using observed relationships between meteorology and stomatal conductance for application to O3. 10 In this study we present a simplified multi-layer atmosphere-biosphere exchange model and investigate the sensitivity of NOx canopy fluxes, ozone production, NOx vertical profiles, and NOx lifetimes to different parameterizations of stomatal conductance and deposition velocity. We consider here both the Wesely model and the similarly simplistic approach of Emberson et al. (2000) that incorporates effects of VPD and SWP. We restrict our considerations to the effects of different stomatal resistance parameterizations on predicted deposition velocities, as the magnitude of the mesophyllic resistance 15 remains uncertain and is assumed to be comparatively small in atmospheric models (Zhang et al., 2002). We also restrict our considerations to NO2 deposition, as NO deposition has been shown to be negligible in comparison (Delaria et al., 2018). There have been many studies investigating the effects of dry-deposition parameterizations on deposition velocities-particularly of ozone-and the abilities of different modeling schemes to reproduce observational data for other molecules such as NO2, NO, H2O2, HNO3, hydroxy nitrates, alkyl nitrates, peroxyacyl nitrates, etc. (Zhang et al., 1996;Wang et al., 1998b;Emberson et 20 al., 2000;Ganzeveld 2002;Buker et al., 2007;Hardacre et al., 2015;Nguyen et al., 2015). However, there has been little evaluation of how changes in dry deposition of NO2 may affect surface mixing ratios and chemistry of important atmospheric species. Assessing the sensitivity to NO2 deposition is crucial not only for evaluating the potential impact of uncertainties in dry-deposition parameterizations for global and regional models, but for understanding how a changing climate may influence NOx, surface ozone, and the nitrogen cycle. 25

Model description
We have constructed a simple atmospheric model for investigating the influence of leaf-level NO2 foliar deposition on canopy scale NOx lifetimes and concentrations. The model consists of eight vertical boxes within the planetary boundary layer (PBL), taken to be 1000 m during the day and 60 m at night . The increase in PBL height during the day is treated as a Gaussian function of time with 98% of the integrated area contained between sunrise and 30 sunset, with the maximum height reached at the time of maximum daily temperature (Fig.1).
In each box, the change in concentration ( ) of species , is calculated using the time-dependent continuity equation: where the terms on the right are the chemical production, chemical loss, emission, deposition, advection, and turbulent flux, respectively. In each box ( =1-8) the altitude (z) is considered as the average of the altitudes at the upper boundaries of boxes and − 1 (the midpoint of box ). The change in concentration for species is calculated for each time step t = 2 s (Table   1). 5 , = ( , + , + , + , + , + , ℎ ) ( 2) where ℎ is the width of box . The only species not treated in this manner is the hydroxyl radical (OH), which is calculated using a steady-state approximation.

Deposition
The deposition flux (F dep ) of each depositing species in the canopy is calculated according to: 10 where LAI is the leaf area index, and is the deposition velocity. The deposition velocities are calculated according to: where R is the total resistance to deposition.
where , , , , and are the aerodynamic, boundary layer, cuticular, stomatal, and mesophilic resistances, respectively. These resistances describe the turbulent transport of a gas to the surface ( ), molecular transport through a thin layer of air above the leaf surface ( ), and deposition to the leaf surface ( ) (Baldocchi et al., 1987). is dependent upon plant physiology and the chemical and physical properties of the deposition compounds. is determined by 20 deposition to the leaf cuticles ( ), diffusion through the stomata ( ), and chemical processing within the mesophyll ( ).
We do not allow for emission of NO or NO2 from leaves, consistent with recent laboratory observations that have observed negligible compensation points for these molecules (Chaparro-Suarez et al., 2011;Breuninger et al., 2013;Delaria et al., 2018).
All boundary, aerodynamic, cuticular, and soil resistances of O3, HNO3, CH2O, alkyl nitrates (ANs) and acylperoxy nitrates (APNs), HC(O)OH, ROOH, and H2O2 are calculated according to . The cuticular and mesophyllic 25 resistances for NO2 and NO are adjustable input parameters. Stomatal resistances are determined from the stomatal conductance to water vapor (gs) calculated using either Eq. 7 (Wesely, 1989), or Eq. 8 (Jarvis et al., 1976;Emberson et al., 2000), hereafter referred to as the Wesely and Emberson schemes, respectively: where is the species-specific maximum stomatal conductance, is a species-specific scaling factor to the minimum stomatal conductance, SR is the solar radiation in W m -2 , and ℎ , , ℎ , , and are functions representing modifications to the stomatal conductance due to leaf phenology, soil water content, irradiance, temperature, and vapor pressure deficit, respectively (Eq 9-12). ) + min (12) T opt and min are the optimal and minimum temperature required for stomatal opening. PPFD is the photosynthetic photon 10 flux density and ℎ is a species-specific light response parameter. min and max are the vapor pressure deficit at which stomatal opening reaches a minimum and maximum, respectively. min and max are the soil water potentials at which stomatal opening reaches a minimum and maximum, respectively. All model calculations represented the peak growing season when ℎ = 1. , , and ℎ were calculated according to Emberson et al. (2000) using parameters found in Table 2. 15

Site description
The model was evaluated with comparison to observations from the Biosphere Effects on Aerosols and  (Table 1) (Fig. 2a). Meteorological conditions and soil NO emissions used in the model simulation were those reported by Min et al. (2014). Diurnal soil water potentials (SWP) were values reported in a geological survey of nearby Sierra sites in a comparatively wet year (Ishikawa and 25 Bledsoe, 2000;Stern et al., 2018).
For UMBS-2012 calculations, the modelled canopy included an overstory height of 20 m with a one-sided LAI of 2.5 m 2 m -2 , and an understory height of 4 m with a LAI of 1 m 2 m -2 (Bryan et al. 2015). Model simulations were run for August 8, 2012 using conditions from the UMBS mixed hardwood forest located in northern Michigan (45°33'32" N, 84°42'52" W) (Table 1) (Fig 2b). Daily temperatures, VPDs, soil NO emissions and site-specific parameters used in the model simulations were those reported in Geddes and Murphy (2014), and Seok et al. (2013).
Temperature and relative humidity used in the model were sinusoidal fits to observations of minimum and maximum daily temperature and relative humidity from the corresponding field measurement site. The relative temperature decrease as a function of altitude was calculated using a fit to observations during BEARPEX-2007, as presented by Wolfe and Thornton 5 (2011). Solar zenith angles (SZA) and photosynthetically active radiation (PAR) were calculated every 0.5 h for each location and time period using the National Center for Atmospheric Research TUV calculator (Madronich and Flocke, 1999) and fit using a smoothed spline interpolation. Within the canopy, extinction of radiation ( ) was calculated following Beer's law: where is the radiation extinction coefficient, is the solar zenith angle,and is the cumulative LAI calculated as the sum of one-half the LAI in box and the total LAI in the boxes above box .

Vertical transport and advection
The turbulent diffusion flux ( ( )) is represented in the model using K-theory, according to the Chemistry of Atmosphere-Forest Exchange (CAFE) Model . 15 where , is the change of concentration in species in box during each timestep and is the difference between the midpoints of boxes and + 1. ( ) above the canopy is based on the values from Gao et al. (1993) and below is a function of friction velocity calculated according to  and is a function of the diffusion timescale ratio ( / )defined as the ratio of the "time since emission" of a theoretical diffusing plume (τ) and the Lagrangian timescale (TL)-and 20 the friction velocity ( * ) . The details of the parameterization of turbulent diffusion fluxes is documented elsewhere  and based on the works of Raupach (1989) and Makar et al. (1999). The height dependent friction velocity (u(z) * ) is attenuated from the above-canopy u * according to Yi et al. (2008). Although Finnigan et al. (2015) identified flaws in this treatment, we believe it is sufficient for our focus on illustrating generalizable qualitative trends. 25 The resulting residence time in the canopy is approximately 2-3 min for model conditions during the day. Our model is a simple parameterization of turbulent processes and as such will only capture mean vertical diffusion. Other works (Collineau and Brunet, 1993a;Raupach et al., 1996;Brunet and Irvine, 2000;Thomas and Foken, 2007;Sörgel et al., 2011;Steiner et al., 2011) have shown that "near-field" effects of individual canopy elements and coherent turbulent structures can play an important role in canopy exchange. These more intricate processes are not captured explicitly by our simple model. 30 Previous work (Gao et al., 1993;Makar et al., 1999;Stroud et al., 2005; have also utilized fairly simple representations of canopy exchange in local and regional models, and as such K-theory is likely sufficient to represent average vertical diffusion for the purposes of our study. Advection in the model is treated as a simple mixing process in each model layer.  Seok et al., 2013) and are used to maintain reasonable background concentrations (Table S1). Concentrations of NOx, O3, and some VOCs at both sites were influenced by emissions from nearby cities and consequently had sources outside the canopy. For the BEARPEX-2009 model runs, the maximum daily advection concentration was reached at around 10 17 hrs, based on field observations of higher NOx plumes from near-by Sacramento in the afternoon Min et al., 2014). The diurnal advection concentrations of NOx were modelled with a sinusoidal function in the range 0.1-0.35 ppb (Table S1). For UMBS all advection concentrations were constant.

Chemistry
Chemistry in the model is based on reaction rate constants from the JPL Chemical Kinetics and Photochemical Data 15 Evaluation No. 18 (Burkholder et al., 2015). Photolysis rates are calculated as a function of solar zenith angle (SZA), which was constructed using a smoothed spline interpolation fit of photolysis rates calculated with the TUV calculator (Madronich and Flocke, 1999) at every ten-degree interval of the zenith angle. The simplified reaction scheme included in the model is based on the model presented in Browne and Cohen (2012). The model includes both daytime and night-time NOx chemistry and a simplified oxidation scheme. In this simplified case, oxidation of volatile organic compounds (VOCs) during the daytime 20 results in the production of peroxy radicals (RO2), treated as a uniform chemical family. To be applicable to a range of forest types, we also include adjustable parameters, kOH and kNO3 for the average rate constant for reaction of VOC with OH and NO3, respectively. kOH and kNO3 are effective values adjusted in the model based on site-specific VOC composition and observations of OH reactivity. A complete list of reactions and rate constants included in the model is shown in Table S2.

BVOC emissions 25
Emissions rates (molecules cm -3 s -1 ) of biogenic volatile organic compounds (BVOCs) in the canopy are calculated via: where (molecules cm(leaf) -2 s -1 ) is the basal emission rate of VOC, ℎ is the total height of the box, and and are corrections for light and temperature (Guenther et al., 1995). 30

Evaluation of NOx fluxes and lifetimes
The model was used to assess the impact of NO2 deposition parameters on the NOx budget, lifetime, loss, and vertical profile within a forested environment. In each box, the rates of NOx loss with respect to nitric acid formation, alkyl nitrate formation, and deposition were calculated from Eq. 17-19.
α is the fraction of the NO + RO2 reaction that forms alkyl nitrates and β is the fraction of the NO3 + BVOC reaction that forms alkyl nitrates. The NOx lifetime was then scaled to the entire boundary layer by summing over the products of the lifetime and boundary layer fraction ( ℎ / ) in each box 10 NOx was treated as the sum of NO, NO2, and all short-lived products, including NO3, 2N2O5, and peroxyacetyl nitrate (PAN) (Romer et al., 2016). Deposition of PAN was not considered.
We also calculated the 24 h average vertical fluxes (Eq. 14) of NOx , and used the flux through the canopy to estimate the fraction of soil emitted NOx ventilated to the troposphere above. Because PAN formed during the nighttime is expected to 15 re-release NOx to the atmosphere during the day, in this calculation, PAN was included as part of the NOx budget.

Sensitivity to parameterizations
We assessed the sensitivity of the model to τ/TL, the radiation extinction coefficient (krad), the aerodynamic leaf width (lw), LAI, soil NO emission (eNO), and α. These parameters are simplifications of complex physical processes and not always easily constrained by observations. The total deposition velocity of NO2 chosen for these assessments was 0.2 cm s -1 during the 20 daytime and 0.02 cm s -1 during the nighttime, based on values of gmax and gmin chosen for Blodgett forest (discussed above) and typical values for deposition velocity observed for a variety of species in the laboratory (Teklemariam and Sparks, 2006;Chaparro Suarez et al., 2011, Breuninger et al., 2013, Delaria et al., 2018. The largest effects were observed for changes in α, LAI, and soil NO emission. LAIos and LAIus were scaled from their values of 1.9 m 2 /m 2 and 3.2 m 2 /m 2 , respectively by a factor of 0.25 and 1.5. Increasing the scaling factor from 0.25 to 25 1.5 resulted in a decrease of NOx lifetimes, above canopy concentration, and average canopy flux of 24%, 27%, and 36%, respectively (Fig. S1). Increasing α from 0.01 to 0.1 resulted in a decrease in NOx lifetimes, above canopy concentrations, and average canopy fluxes of 75%, 38%, and 39%, respectively (Fig. S2). For all other model runs an α of 0.075 was chosen, in accordance with observations from regions primarily influenced by BVOCs (eg. monoterpenes, isoprene, 2-methyl-3-buten-2-ol). Increasing the maximum soil NO emission from 1 to 10 ppt m s -1 increased the in-canopy enhancement from 28% to 30 140% relative to above-canopy NOx concentrations (Fig. S3b). The fraction of soil-emitted NOx ventilated through the canopy also increased from 45% to 64% (Fig. S3a). The large effect of soil NO emission on NOx fluxes implies that this highly variable parameter (Vinken et al., 2014) is also important to constrain in chemical transport models. Further discussion of soil NO emission is, however, beyond the scope of this study.
Very small effects on NOx were observed for changes in the parameters τ/TL, krad, or lw. The minor changes caused by variations in these parameters are listed below for completeness: 5 τ/TL represents the diffusion timescale ratio, a full description of which can be found in .
Larger τ/TL represents faster diffusion and vertical transport within the canopy layer, and shorter residence times in the canopy.
We find that altering this parameter from 1.2 to 8 (representing a change in residence time from 650 s to 62 s) caused a 9.9%, 4.4%, and 8.7% increase in average canopy fluxes, NOx lifetimes and above canopy concentration, respectively (Fig. S4). For all subsequent model runs, a value of 2 for τ/TL was chosen, resulting in a canopy residence time during the day of 152 s and 10 194s for Blodgett Forest and UMBS, respectively, calculated using Eq.21.
The boundary layer resistance, or laminar sublayer resistance, Rb, is dependent upon the aerodynamic leaf width, lw where ν = 0.146 cm 2 s -1 is the kinematic viscosity of air, D is the species-dependent molecular diffusion coefficient, c is a tunable constant set to 1 for this study, and u * (z) is the height-dependent friction velocity that is a function of u* and LAIcum . lw depends upon the vegetation species. A value of 1 cm was chosen for the overstory and 2 cm for the understory, as these widths are characteristic of pine trees and understory shrubs in a poderosa pine forest . Species with rapid deposition to the cuticles or the stomata are expected to be more sensitive to errors in lw, 20 such as HNO3 or H2O2. An increase in NOx lifetime, average canopy flux, and above canopy concentration of 1.4%, 2.4%, and 2.8%, respectively, was predicted for a change in lw scaling factor from 0.1 to 2 (Fig. S5). These changes are expected to be greater in forests with a larger average deposition velocity, where Rb makes a greater contribution to the total resistance.
The rates of stomatal gas exchange and photolysis are regulated by the intensity of light that penetrates the canopy.
The extinction of radiation by the canopy, treated as a Beer's Law parameterization (Eq. 113) is exponentially proportional to 25 the radiation extinction coefficient, krad. krad ranging from 0.4-0.65 has been measured for coniferous forests and understory shrubs . The NOx lifetime increased by 2.7% and the canopy fluxes, and above-canopy concentrations decrease by 0.7% and 0.6%, respectively, for a change in krad from 0 to 0.6 (Fig. S6). This effect is expected to be greater for forests with larger LAI. The minimal effect of krad on model results was also observed for multiple canopy profile shapes of equivalent LAI. 30

Model validation: comparison to field observations
To evaluate the applicability of our 1D multilayer canopy model for predicting NOx concentrations and vertical fluxes in a variety of forest environments, we compared the model to observations from BEARPEX-2009 and UMBS-2012. Parameters used in each calculation are shown in Table 1. The model was run using both the Emberson and Wesely stomatal conductance 5 models. Parameters for temperature, drought stress, and maximum and minimum stomatal conductances used in the Emberson model were input for the dominant tree species in the region (Table 2). At the BEARPEX-2009 site, the dominant tree species was ponderosa pine. For this site, and parameters for and were inferred from ponderosa pine stomatal conductance data (Kelliher et al., 1995;Ryan et al., 2000;Hubbard et al., 2001;Johnson et al., 2009;Anderegg et al., 2017), and ℎ was inferred from measurements of the canopy conductance during BEARPEX-2009 (Fig 3a). was 10 represented by observations for Scots pine (Altimir et al., 2004;Emberson et al., 1997;Buker et al., 2012) and validated with comparison to stomatal conductance measured via sap-flow during BEARPEX-2009 (Fig 3a). At UMBS the dominant species are quaking aspen and bigtooth aspen, with many birch, beech, and maple species also present (Seok et al., 2013). Data for a European beech tree species was used to represent stomatal conductance parameters (Buker et al., 2007;Buker et al., 2012) and SWP stress (Emberson et al., 2000). These parameters were validated with comparison to stomatal conductance calculated 15 from water vapor and latent heat flux measurements during UMBS-2012 using an energy-balance method according to Mallick et al. (2013) (Fig 4a).
The model replicates key features of the canopy fluxes and above-canopy NOx daytime mixing ratios from the 2009 BEARPEX campaign (Fig.3). The average daytime above-canopy NOx mixing ratios during the duration of BEARPEX-2009 was 253 ppt, with observations ranging from 80-550 ppt of NO2 and 10-100 ppt of NO (Min et al., 2014). The general daily 20 trends in observations of NOx mixing ratios are captured by both the Wesely and Emberson cases-with minimum NOx mixing ratios occurring in the late morning, an increase of NOx in the afternoon, and maximum NOx concentrations of 450-500 ppt reached in the evenings, primarily as a result of high-NOx plumes from near-by Sacramento in the afternoon Min et al., 2014) (Fig. 3b). However, both model scenarios predict a slower than observed decrease in NOx mixing ratios from the evening to the early morning, larger mid-morning fluxes than observed (by ~0.5-1.5 ppt m s -1 ), and fail to represent 25 the in-canopy enhancement of NOx (~50 ppt), relative to above-canopy mixing ratios, observed in the evening (Fig 3). The above-canopy vertical NOx flux predicted in both model cases also agrees reasonably well with observations, with the Emberson case representing morning and midday NOx fluxes slightly better than the Wesely case. This relatively good agreement between the Emberson case and observed fluxes is also demonstrated in Fig 3d by the agreement between modelled and observed canopy NOx enhancements. There is, however, generally little difference between Emberson and Wesely model 30 cases for this site during the period considered (Fig 3). This is likely due to the good agreement in both the Emberson and Wesely cases to observations of stomatal conductance (Fig 3a).
We also observe similar correspondence between the model and key features of the UMBS-2012 observations (Fig   4). NO and NO2 mixing ratios and canopy fluxes are both within the range of observations. The model predicts a maximum of ~40% lower NO2 in the morning and ~30% higher NO2 at night than what was observed (Fig 4b). Differences between the Wesely model and Emberson model were negligible for this site. This is likely due to a higher humidity in the summer in this region and larger soil moisture, reducing the prediction for midday and late afternoon VPD stress by the Emberson model, as 5 can be seen by the similarity in the predicted gs by the Emberson and Wesely models (Fig 4a).

Effects of maximum stomatal conductance
The BEARPEX-2009 case was simulated using the Wesely model for different values of the maximum stomatal conductance ( ) (Fig 5), with advection concentrations of NOx set to zero. The range of currently represented in the literature during peak growing season for forested regions ranges from 0.2-0.8 cm s -1 (Kelliher et al., 1995;Emberson et al., 10 1997;Emberson et al., 2000;Ryan et al., 2000;Hubbard et al., 2001;Altimir et al., 2003;Fares et al., 2013). This range reflects differences in forest types and a wide variety of tree species. Global CTMs using the Wesely parameterization currently include of 1.4, 0.77, and 1 cm s -1 for deciduous, coniferous, and mixed forests, respectively (Wesely, 1989;Wang et al., 1998a). NOx is ventilated through the canopy with no foliar deposition ( = 0 cm s -1 ). In contrast, 44% of soil-emitted NOx is taken 15 up by the forest and 56% ventilated through the canopy when the maximum stomatal conductance is 1.4 cm s -1 . Figures 5c and   5d show the effects of on the diurnal flux through the canopy and the diurnal above canopy NOx mixing ratio, respectively. Compared with no foliar deposition, a of 1.4 cm s -1 results in ~60% reduction in the canopy flux and ~50% reduction in the above-canopy NOx mixing ratio at noon. (Fig. 5c, d). In Figure 6a we show the fraction of soil-emitted NOx ventilated through the canopy as a function of . The model suggests a maximum foliar reduction of NOx of ~60% for a 20 canopy of 10 m and total LAI of 5.1 m 2 /m 2 . Our model also predicts that changes in have a greater overall impact on canopy NOx fluxes at larger leaf resistances and slower foliar uptake. In the range for of ~0-0.5 cm s -1 , variation in can have a large impact on the predicted canopy fluxes of NOx, which would in turn have a large impact on concentrations and fluxes of O3. These values of results in deposition velocities in the range expected for most forests, based on laboratory measurements of leaf-level deposition (Hanson and Lindberg, 1991;Rondon and Granat, 1994;Hereid and Monson, 2001;25 Teklemariam and Sparks, 2006;Pape et al., 2008;Chaparro-Suarez et al., 2011;Breuninger et al., 2013;Delaria et al., 2018) and global analysis suggesting 20-50% reductions in soil-emitted NOx by vegetation (Jacob and Wofsy, 1990;Yienger and Levy, 1995). Model calculations also predict a strong effect on the lifetimes of NOx, as shown in Figure 6b, with maximum stomatal conductances of 0.1 cm s -1 , 0.3 cm s -1 , and 1.4 cm s -1 reducing the NOx lifetime by ~ 0.7 hrs (~7%), ~1.8 hrs (~18%), and ~3.6 hrs (~36%), respectively compared with no deposition. Similar trends (not shown) were also observed using 30 parameters for UMBS.

Emberson model vs. Wesely model comparison
The relative importance of including parameterizations of VPD and SWP in the calculation of stomatal conductance and overall deposition velocity is expected to be regionally variable, along with regional variations in dominant tree species, soil types, and meteorology. We ran the model using BEARPEX-2009 conditions using both the Wesely and Emberson stomatal conductance models under "dry" and "wet" conditions. Under the "dry" scenario the SWP daily minimum and maximum were 5 -2.0 MPa and -1.7 MPa, respectively, with the daily minimum reached at sunset. A minimum daily RH of 40% occurred at noon, with a maximum at midnight of 65%. Summertime is often even drier in regions of the western United States, so these "dry" parameters are conservative estimates for many forests. Under the "wet" scenario the SWP daily minimum and maximum were -0.5 MPa and -0.1 MPa, respectively. The maximum and minimum RH were 90% and 80%, respectively. The values for soil moisture and relative humidity chosen were based on observations of SWP by Ishikawa and Bledose (2000) and the long-10 term climate data record at Auburn Municipal Airport (38.9547° N, -121.0819° W) from NOAA National Centers for Environmental Information.
The results of the Wesely and Emberson "wet" and "dry" model runs are shown in Figure 8. There was only a slight decrease of the in-canopy NOx enhancement and of the canopy fluxes in the Wesely "wet" case, presumably due to a slight increase in OH radicals at higher RH. Predictably, the difference in the modelled deposition velocities was quite dramatic 15 between the Emberson "wet" and "dry" cases. In the "dry" scenario, the deposition velocity reached a maximum in the late morning, but rapidly declined after noon. The maximum deposition velocity reached was also substantially reduced (Fig 7a).
Using the "wet" Emberson stomatal conductance model, the NOx flux out of the forest was reduced by 16% midday compared to the "dry" case, and the percent of soil NOx removed within the canopy was increased from 18% to 30% (Fig 7). The model calculates a substantial impact on above-canopy NOx mixing ratios (Fig. 8), with a maximum of ~30% difference in NOx in 20 the afternoon between "wet" and "dry" days using the Emberson parameterizations, compared with ~10% difference using the Wesely model. Using the Emberson parameterization of stomatal conductance, deposition during "wet" days is predicted to contribute substantially more to the total NOx loss (~40%), with only ~15% contribution predicted for "dry" days ( Fig. 9). Under the Wesely model, where stomatal conductance is parameterized only with temperature and solar radiation, the predicted deposition velocity would be nearly identical between the spring and fall in the western United States and similar 25 semi-arid regions (with comparatively minor temperature effects). While the Emberson model predicts large seasonal differences, the Wesely model fails to account for the dramatic decrease in stomatal conductance seen in the dry seasons in such regions caused by significant reductions in relative humidity and soil water potential (Prior et al., 1997;Panek and Goldstein, 2001;Chaves, 2002;Beedlow et al., 2013). We recognize that the multibox model presented in this work is a simplified representation of physical processes, and as such is not likely to (and is not intended to) provide quantitative 30 exactitude for the trends described above. However, we argue for the necessity of incorporating these conceptual advances for accurately representing canopy processes and predicting their effect on the NOx cycle.

Implications for modelling NO2 dry deposition
As in our multilayer canopy model, the most common current method of parameterizing stomatal and cuticular deposition in large-scale chemical transport models (CTMs) is through the resistance model framework of Baldocchi (1987).
Many global (e.g. WRF-Chem and GEOS-Chem) and regional (e. g. MOZART and CAMx) CTMs calculate the stomatal 5 component of the total deposition resistance using the representation of Wesely (1989), where stomatal conductance is dependent only on the type of vegetation, temperature, and solar radiation. The limitations of this parameterization have been highlighted by observations of a strong dependence of foliar deposition on soil moisture and vapor pressure deficit (VPD) (Kavassalis and Murphy, 2017;Rydsaa et al., 2016). Inadequate descriptions of vegetative species, soil moisture, drought stress, etc., can have a dramatic impact on model results, and result in significant discrepancies between models and 10 observations (Wesely and Hicks, 2000). Failure to account for effects of plant physiology on deposition may result in misrepresentation of deposition velocities, which, as we demonstrate, can have a substantial impact on NOx lifetimes and mixing ratios above and within a forest canopy. This effect will be especially pronounced in areas, such as much of the western United States, where there are frequent periods of prolonged drought. Parameterizations of stomatal conductance, such as those presented in Emberson et al. (2000) and incorporated into some regional-scale CTMs (e.g. EMEP, MSC-W, and CHIMERE), 15 if incorporated into global atmospheric models, could more accurately reflect the dependence of foliar deposition on meteorology and soil conditions. However, additional laboratory and field measurements on diverse plant species are also needed to determine appropriate, ecosystem-specific inputs to these parameterizations.
It should be noted that there have been significant recent advances in optimization approaches of stomatal modelling based on the theory that stomata maximize CO2 assimilation per molecule of water vapor lost via transpiration (Medlyn et al., 20 2011;Bonan et al., 2014;Franks et al., 2017;Miner et al., 2017;Franks et al., 2018). Medlyn et al. (2011) reconciled the empirical widely utilized Ball-Berry model with a theoretical framework optimizing ribulose 1,5 bisphosphate (RuBP) regeneration-limited photosynthesis. However, such methods of water use efficiency optimization do not account for stomatal closure as a result of soil moisture stress. Bonan et al. (2014) further developed a model considering water use efficiency optimization and water transport between the soil, plant, and atmosphere. Such parameterizations are utilized in the 25 Community Land Model (CLM)-a land surface model often incorporated into regional and global climate-chemistry models (Lombardozzi et al., 2015;Kennedy et al., 2019). Although this model provides a physiological and mechanistic basis for stomatal behaviour, it is heavily parameterized, relying on inputs of plant and soil parameters that could be expected to vary significantly across ecosystem types. For this reason, we view these methods as aspirational for incorporation into atmospheric global CTMs. We find the relative simplicity of the Emberson approach more useful for the purpose and scope of parameters 30 for large-scale atmospheric models.

Implications for modelling ozone
NO2, as well as O3, deposition budgets are frequently calculated through inferential methods whereby the deposition velocity is constrained with ambient observations (Holland et al., 2005;Geddes and Murphy, 2014). These inferential models are often complicated by the fast reaction of the NO2-NO-O3 triad, making it difficult to separate chemical and physical processes. Further, these inferential models for determining dry deposition constrained with observations of chemical 5 concentrations and eddy covariance measurements of fluxes are difficult to interpret because of similar chemical and turbulent timescales (Min et al., 2014;Geddes and Murphy, 2014). Emission of NO from soils, rapid chemical conversion to NO2, and subsequent in-air reactions of NOx must be evaluated accurately in in order to correctly infer NOx and O3 atmosphere-biosphere exchange from observations. Our multilayer canopy model applies a simple method of representing these processes and evaluating the separate effects of chemistry and dry deposition on the NOx budget in forests. 10 Since the foliar deposition of NO2 reduces the NOx lifetime and NOx that is transported out of the canopy, it will also reduce the amount of ozone that is produced both within and above the canopy. Ozone production efficiency (OPE) in the canopy is calculated using Eq.23-25: where (O 3 ) is the ozone production rate and (NO x ) is the NOx loss rate. The effect of stomatal conductance to NO2 on OPE is shown in Figure 6c. An increase in from 0 to 0.3 cm s -1 results in a decrease in OPE for the PBL from 24.0 to 20.7 (~14%), and a decrease to 17.0 (~30%) if is 1.4 cm s -1 . This is similar to OPE calculations that have been reported for forests and environments with NOx mixing ratios less than 1 ppb and heavily influenced by BVOC emissions (Marion et al., 20 2001;Browne and Cohen, 2012;Ninneman et al., 2017).
NO2 deposition and the in-canopy chemistry of NO2-NO-O3 also impacts O3 production and removal. O3 deposition is frequently inferred from measurements of O3 concentrations or eddy-covariance measurements (Wesely and Hicks, 2000;Kavassalis and Murphy, 2017). However, because NO2 has a direct impact on ozone production, deposition of NO2 can affect inferences of O3 deposition from observations. The 14% reduction of OPE and more than a 20% reduction in daytime NOx 25 resulting from an increase in from 0 to 0.3 cm s -1 can cause a parallel decrease in O3 concentrations and fluxes independent from O3 chemical loss or deposition. Thus, deposition of NO2 must be taken into account when evaluating O3 deposition losses from observed canopy fluxes.

Implications for near-urban forests
The analysis above suggests that the relative importance of chemical sinks and deposition will vary with NOx concentration. 30 To evaluate the relative importance of NO2 foliar deposition and chemistry as a function of NOx mixing ratio, a simplified 1box model was also constructed with a simplified reaction scheme (Table S3), VOC reactivity of 8 s -1 , of 0.075, and a HOx (HOx ≡ OH + HO2) production rate ( ) of 2×10 6 molecules cm -3 s -1 (similar to conditions observed at BEARPEX-09). RO2, OH, and HO2 were solved for steady-state concentrations and NOx loss pathways were calculated via Eq. 26-29.
where ℎ is the canopy height (15m), ℎ is the planetary boundary layer height (1000 m), and LAI is 5 m 2 /m 2 .
The results from this simplified box model are shown in Figure 9 and agree well with our 1D multi-box model near 10 ppb NOx ( Fig S7). With deposition set to zero, nitric acid formation becomes a more significant sink of NOx than alkyl nitrate 10 formation at around 1 ppb, and nitric acid formation accounts for greater than 70% of the total loss at 100 ppb. With a deposition pathway included, deposition acts as the dominant NOx sink above 30 ppb and at 10 ppb deposition and AN formation are each 20% of the NOx sink. Deposition is approximately 10% of the sink over a wide range of concentrations. Forests in close proximity to urban centers ( Fig S9) may result in a substantial local decrease in NOx (Fig S8, Fig 10). Although the influence of urban or near-urban trees on NOx concentrations would be heavily dependent on meteorological factors (i.e. wind speed and 15 direction), proximity to emission sources, and LAI, it may have some importance on a local or neighborhood scale. This effect may be relevant for understanding and predicting the effects of NOx reduction policies within and near cities. It may also be useful in considering as a direct nitrogen input to the biosphere, not mediated by soil processes.

Conclusions
We have constructed a 1D multi-box model with representations of chemistry and vertical transport to evaluate the impact of 20 leaf-level processes on canopy-scale concentrations, lifetimes, and canopy fluxes of NOx. Our model is able to closely replicate canopy fluxes and above-canopy NOx daytime mixing ratios during two field campaigns that took place in a Sierra Nevada pine forest (BEARPEX-2009) and a northern Michigan mixed hardwood forest (UMBS-2012). We conclude that the widely used canopy reduction factor approach to describing soil NOx removal from the atmosphere within plant canopies is consistent with a process-based model that utilizes stomatal uptake and we recommend that the CRF parameter be replaced with stomatal 25 models for NO2 uptake.
We demonstrate with our 1D multi-box model that NO2 deposition provides a mechanistic explanation behind canopy reduction factors (CRFs) that are widely used in CTMs. We predict a maximum of ~60% reduction in the fraction of soilemitted NOx ventilated through the canopy when stomatal conductances are greater than 0.075 cm s -1 , consistent with the range of global CRFs used in current CTMs (Jacob and Wofsy, 1990;Yienger and Levy, 1995). Our model also predicts that changes 30 in have a greater overall impact on canopy NOx fluxes at larger leaf resistances to uptake (slower foliar uptake). In the range for of ~0-0.5 cm s -1 , errors or variability in stomatal conductance can have a large impact on the predicted canopy concentrations and fluxes of NOx, which would in turn have large impact on concentrations and fluxes of O3. This range of deposition velocities describes the range of uptake rates measured for many tree species and forest ecosystems (Hanson and Lindberg, 1991;Rondon and Granat, 1994;Hereid and Monson, 2001;Teklemariam and Sparks, 2006;Pape et al., 2008;Chaparro-Suarez et al., 2011;Delaria et al., 2018). Model calculations also predict a similar trend on the lifetimes of NOx, with 5 a maximum reduction in the NOx lifetime by ~4 hrs (>40%) compared with no deposition.
The large effect that small changes in stomatal conductance can have on NOx lifetimes, concentrations, budget, and O3 production makes it very important to accurately parameterize in atmospheric models. Most global scale chemical transport models parameterize stomatal conductance using the representation developed by Wesely (1989) (Jacob and Wofsy, 1990;Wang and Leuning, 1998;Ganzeveld et al., 2002;Verbeke et al., 2015). These do no account for the effects of VPD, SWP, 10 CO2 mixing ratio, or other factors known to influence stomatal conductance (Hardacre et al., 2015). We show that incorporating vapor pressure deficit and soil water potential-using the parameterization of Emberson et al. (2000)-has a substantial impact on predicted NO2 deposition, with the percent of soil NOx removed within the canopy increasing from 18% to 30% in wet (low VPD and high SWP) conditions compared to dry conditions in the location of BEARPEX-2009. Under the Wesely model, where stomatal conductance is parameterized only with temperature and solar radiation, the predicted deposition velocity 15 would be nearly identical between "wet" and "dry" days and between the spring and fall in semi-arid regions (e.g. much of the western United States, the Mediterranean Basin, the west coast of South America, parts of northwest Africa, parts of western and southern Australia, and parts of South Africa). The dominant effect of stomatal opening on NO2 deposition causes an important time of day and seasonal behaviour that should be extensively explored with observations of NOx fluxes and concurrent models to confirm the role of deposition in a wider range of environs and more thoroughly vet the conceptual model 20 proposed here.
Code availability. The model reported in this paper has been deposited in GitHub, https://github.com/erd02011/NOxmodel_ACP2019. 25 Data availability. The data presented in this study from Blodgett Forest and the University of Michigan Biological field Station have been published in earlier work by Min et al. (2014) and Geddes and Murphy (2014), respectively.
Author contributions. ERD built the model, preformed data analysis, and prepared all figures. ERD wrote the manuscript in consultation with RCC. RCC supervised the project. 30 Competing interests. The authors declare that they have no conflict of interest.
University for providing data from the UMBS field site; and J. Geddes for constructive comments that improved the manuscript.   Ishikawa and Bledsoe (2000) and Stern et al. (2018)