Evaluation of and updates to the oxidized reactive nitrogen gaseous dry-deposition parameterization from the GEOS-Chem model, including a pathway for ground surface NO₂ hydrolysis

Dry deposition is a major loss pathway for reactive nitrogen species from the atmospheric boundary layer. We evaluate isolated components of the parameterization for species-specific gaseous dry-deposition velocity V_d(x) for HNO₃ and NO₂ from the GEOS-Chem chemical transport model by running a stand-alone version of V_d code in single-point mode to enable a more direct comparison to field observations. Improved measurement–model agreement results mainly from (i) updates to the calculation of molecular diffusivities and (ii) the representation of ground surface NO₂ hydrolysis in the formulation of non-stomatal uptake. We evaluate the parameterization for non-stomatal dry deposition of NO₂ by comparing to eddy-covariance-inferred nocturnal V_d(NO₂) over Harvard Forest. We address a large low bias (−80 %) in simulated nocturnal V_d(NO₂) by representing NO₂ heterogeneous hydrolysis on deposition surfaces, paying attention to chemical flux divergence, soil NO emission, and canopy surface area effects. Finally, we evaluate the updated oxidized reactive nitrogen (NO_y) dry-deposition parameterization by comparing to eddy-covariance-inferred V_d(NO_y) over Harvard Forest, finding that a modest nocturnal low bias (−19 %) remains in simulated V_d(NO_y) due to the compensating effects of updates to the calculation of molecular diffusivities (28 % reduction in nocturnal V_d(NO_y)) and the representation of NO₂ heterogeneous hydrolysis (25 % increase in nocturnal V_d(NO_y)). These developments are a first step towards a tractable representation of NO₂ hydrolysis in a dry-deposition scheme and have important implications for the near-surface NO₂ lifetime through a mechanism involving HONO emission.

1 Introduction

Chemical species constituting oxidized reactive nitrogen (NO_y) form a major component of atmospheric reactive nitrogen (N_r≡ NO_y+ reduced nitrogen species), which together play a central role in atmospheric chemistry by modulating the oxidative capacity of the atmosphere through nitrogen oxides (NO_x≡ NO + NO₂) (Crutzen, 1979), contributing to nitrogen loading of natural ecosystems (Clark et al., 2018) and influencing air quality (Fields, 2004). Accurate knowledge of the sources and sinks of N_r is vital for understanding and modeling atmospheric chemistry, including the sensitivity of air quality to changes in anthropogenic emissions. Dry deposition of N_r from the atmospheric boundary layer is an important removal process, typically contributing between one-third and two-thirds of total (wet + dry) deposition (Flechard et al., 2011; Hanson and Linderg, 1991; Munger et al., 1998; Sparks et al., 2008; Walker et al., 2020), but questions remain about its representation in chemical transport models (CTMs).

The atmosphere–surface exchange of N_r may be measured directly via micrometeorological techniques (Businger, 1986; Walker et al., 2020) or under more controlled conditions via enclosure techniques (Breuninger et al., 2012; Hanson and Linderg, 1991). Direct measurements of the above-canopy air–surface exchange of N_r, including via the eddy covariance technique, are technically complex and resource intensive, resulting in a scarcity of flux observations across representative land types and seasons (Walker et al., 2020). Therefore, studies of above-canopy dry deposition tend to be intensive in nature and are typically designed to characterize exchange processes rather than to monitor long-term deposition patterns. Dry-deposition budgets thus fall in the realm of inferential methods, where deposition fluxes F_x are inferred from parameterizations of above-canopy deposition velocity V_d – a first-order rate coefficient for heterogeneous surface reaction/uptake for a specific gas x to a specific bulk surface/land type from a specified height:

$$\begin{array}{}\text{(1)}& F_x=-V_{\mathrm{d}}\left(x\right)\left[x\right].\end{array}$$

By convention, downward fluxes toward the surface are negative values represented by positive deposition velocities. N_r component concentrations [x] from which dry-deposition budgets may be inferred have been obtained from (i) surface networks such as the US CASTNET (Clarke et al., 1997) and Canadian CAPMoN (Zhang et al., 2009), (ii) chemical transport models (Dennis et al., 2013; Zhang et al., 2012, 2018), and (iii) satellite observations (Geddes and Martin, 2017; Kharol et al., 2018; Nowlan et al., 2014).

Deposition velocity represents a bulk quantity with contributions from complex processes including turbulent and molecular diffusion in air; meteorological influence on the physical, chemical, and biological state of surfaces; and species-specific interfacial chemistry. The most common parameterization of V_d in large-scale CTMs considers the deposition pathway as a series of three resistances (Baldocchi et al., 1987; Wesely and Hicks, 1977):

$$\begin{array}{}\text{(2)}& V_{\mathrm{d}}(z,x)={\displaystyle \frac{\mathrm{1}}{R_{\mathrm{a}}\left(z\right)+R_{\mathrm{b}}\left(x\right)+R_{\mathrm{c}}\left(x\right)}},\end{array}$$

where, for bulk-canopy V_d above a projected ground area, R_a(z) is the aerodynamic resistance to turbulent transport from a specified height z and is common for all species, R_b(x) is the species-specific quasi-laminar boundary layer resistance to transport through the thin non-turbulent layer in direct contact with surfaces, and R_c(x) is the bulk-canopy surface resistance for a specific land type. Expressions for R_a and R_b can be obtained from micrometeorological flux–gradient relationships (Garratt, 1992; Wesely and Hicks, 1977) and vary as a function of surface roughness, wind speed, diabatic stability, and molecular diffusivity in air. For highly soluble species such as HNO₃ and H₂O₂, contributions from R_c are nominally small, with the resulting deposition varying between being R_a and R_b limited depending on the state of the turbulence (Nguyen et al., 2015). For species with low aqueous solubility or limited interfacial reactivity, R_c is the limiting term, except under very stable conditions (Toyota et al., 2016). Given the complexity and variability in canopy types and species-specific surface reactivities, R_c is difficult to treat theoretically, with parameterizations relying heavily on empirical formulations.

The most common parameterization of R_c used by large-scale atmospheric models, including the widely utilized WRF-Chem and GEOS-Chem CTMs, is the Wesely 1989 algorithm (Wesely, 1989; hereafter referred to as W89) or modifications thereof (Hardacre et al., 2015). In this scheme, the bulk canopy is treated as a single uniform surface or “big leaf”, with stomatal and various non-stomatal deposition pathways acting in parallel. Trace-gas-specific component surface resistances are calculated following basic similarity relations, including solubility relative to SO₂ and oxidative potential relative to O₃. Zhang et al. (2003a) present a parameterization of R_c(x) for use in air quality models, including at the global scale, using similarity arguments to SO₂ and O₃ as was done in W89, with updates including online computation of within-canopy aerodynamic resistance, the influence of leaf water vapor pressure deficit and water stress on stomatal resistance, and updated parameterizations of non-stomatal surface resistances for O₃ (Zhang et al., 2002b) and SO₂ (Zhang et al., 2003b). Zhang et al. (2003a) note that application of the algorithm (hereafter referred to as Z03) to compounds for which few to no deposition flux observations exist will continue to be a source of uncertainty. Bulk-canopy surface resistances deviating from W89 similarity to SO₂ and O₃ have been observed for NO₂ (Eugster and Hesterberg, 1996; Horii et al., 2004; Stocker et al., 1995), PAN (Shepson et al., 1992; Sun et al., 2016; Turnipseed et al., 2006), and many other species (Nguyen et al., 2015). Wu et al. (2012) compare observed V_d(PAN) over a coniferous forest to deposition velocities parameterized according to both the W89 (WRF-Chem) and Z03 (NOAH) schemes and find underestimates greater than a factor of 2, motivating efforts to fit non-stomatal R_c(PAN) directly from above-canopy flux observations. Using the eddy covariance flux dataset from Nguyen et al. (2015), the Z03 scheme was extended by Wu et al. (2021) to additional species by fitting non-stomatal uptake of oxidized volatile organic compounds (VOCs) and hydrogen cyanide directly from observations. Wu et al. (2021) maintain the Z03 algorithm structure through similarity to SO₂ and O₃; however, they suggest that future developments to dry-deposition schemes should consider other species-specific processes and reactions affecting measured uptake, including below-sensor chemical flux divergence, enzymatic reactions, and other non-stomatal processes/reactions.

A main result of the Horii et al. (2004) analysis of an extensive eddy covariance flux dataset of NO₂ over a northeastern US mixed forest (Harvard Forest) from April to November was that a persistent deposition process was active at night, yielding NO₂ deposition velocities on average of ∼ 0.2 cm s⁻¹, with values up to 0.5 cm s⁻¹ noted under high-NO₂ loads of ∼ 30 ppb. This observation is contrary to the widely used W89 parameterization, which does not allow significant surface uptake of NO₂ at night when leaf stomata are assumed to be closed or during vegetatively dormant seasons. Geddes and Murphy (2014) monitored eddy covariance fluxes of NO, NO₂, and NO_y above midlatitude (∼ 45° N) summertime mixed hardwood forests in Ontario (Canada) and Michigan (USA), finding that, on average, NO_x fluxes were indistinguishable from zero for these relatively low-NO_y environments (< 2 ppb on average). However, infrequent nocturnal events with high NO_x $/$ NO_y ratios and large downward NO_y fluxes could be interpreted as yielding NO₂ deposition velocities similar to the average values of Horii et al. (2004). Geddes and Murphy (2014) were careful to note that above-canopy fluxes of NO_x are influenced by not only deposition processes but also within-canopy emissions and chemistry, resulting in above-canopy fluxes of NO_x that are confounded by a combination of counteracting mechanisms, which render flux observations difficult to interpret. Horii et al. (2004) considered below-sensor chemical flux divergence of NO₂ due to the formation and subsequent hydrolysis of N₂O₅, where the maximum rate of loss was insufficient to account for the observed downward nocturnal NO₂ flux. To reconcile this, they proposed a non-stomatal hydrolysis pathway for uptake of NO₂ on ground and canopy surfaces – a reaction that has been suspected to be of atmospheric relevance in the field for some time (Harrison and Kitto, 1994; Harrison et al., 1996).

The hydrolysis of NO₂ on hydrated surfaces is a well-known heterogeneous reaction from lab investigations, yielding adsorbed HNO₃ and evolved nitrous acid (HONO):

$$\begin{array}{}\text{(R1)}& \mathrm{2}{\mathrm{NO}}_{\mathrm{2}\left(g\right)}+{\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\left(\mathrm{ads}\right)}\stackrel{\mathrm{surface}}{⟶}{\mathrm{HONO}}_{\left(g\right)}+{\mathrm{HNO}}_{\mathrm{3}}\left(\mathrm{ads}\right).\end{array}$$

Despite the stoichiometry of Reaction (R1), first-order kinetics for NO₂ have generally been observed in the lab and field, with a rate dependence on surface area density (as expected for collision-limited heterogeneous catalysis), surface water content, and other surface chemical properties (Finlayson-Pitts et al., 2003; Finlayson-Pitts, 2009; Lammel, 1999; Spataro and Ianniello, 2014). In addition to hydrated ground (Kurtenbach et al., 2001; Lammel, 1999; Ren et al., 2020; VandenBoer et al., 2013) and aerosol (Bröske et al., 2003; Burkholder et al., 2015; Crowley et al., 2010; Tan et al., 2016) surfaces, Reaction (R1) has been suggested to be occuring on the sea surface (Wojtal et al., 2011; Yang et al., 2021; Zha et al., 2014), on snow and ice surfaces (Beine et al., 2001; Kim and Kang, 2010), and on indoor surfaces (Collins et al., 2018; Febo and Perrino, 1991; Spicer et al., 1993). Spicer et al. (1993) and Collins et al. (2018) both found an indoor lifetime of NO₂ to reactive loss (HONO producing) on residential interior surfaces on the order of 1 h in well-mixed air – lower than typical ambient NO₂ chemical lifetimes, which are on the order of hours in regional (Kenagy et al., 2018; Shah et al., 2020) or urban outflows (Laughner and Cohen, 2019) and remote-forest environments (Browne and Cohen, 2012). Reaction (R1) may be an especially important surface removal process during summertime nights or winter months when NO₂ is longer lived, with lifetimes on the order of 10 h to more than a day (Browne and Cohen, 2012; Kenagy et al., 2018; Martin et al., 2003). Reaction (R1) has also been implicated in the uptake of NO₂ through leaf stomata, where it may be an important contributor to NO₂ deposition within the moist and high-surface-area substomatal cavities (apoplast) of leaves (Ammann et al., 1995).

Despite the evidence for Reaction (R1) occurring on nearly any environmental surface with adsorbed water, regional and global CTMs have yet, to our knowledge, to update dry-deposition parameterizations of NO₂ to include this effect, potentially underestimating and/or misrepresenting V_d(NO₂) at night and throughout vegetatively senescent periods when stomatal uptake would be weak or absent. In this study, we compare simulated dry-deposition velocities from the GEOS-Chem CTM to above-canopy observations of V_d(NO₂) and V_d(NO_y) inferred from an extensive publicly available dataset of NO₂ and NO_y eddy covariance fluxes and speciated NO_y concentration measurements over Harvard Forest (Munger and Wofsy, 2023), paying attention to soil NO emission, chemical flux divergence, and canopy surface area effects. Prior to updating simulated R_c(NO₂) to include Reaction (R1), we conduct sensitivity tests to evaluate the parameterization of R_a and R_b used in GEOS-Chem by comparing to daytime deposition velocities of rapidly depositing species inferred by the method of eddy covariance over a southern US temperate forest (Nguyen et al., 2015). Specifically, we comment on the effects that site-specific roughness length, reference height, and the roughness sublayer have on the simulation of daytime R_a, followed by the correction of a positive bias in calculated molecular diffusivities that greatly improves the simulation of daytime V_d(HNO₃) via a large relative increase in R_b.

2 Reference model and measurements

2.1 Reference algorithms for computing gaseous dry-deposition velocities

2.1.1 GEOS-Chem dry-deposition module

To facilitate site-specific comparisons to measured deposition velocities, we use a stand-alone version of the gaseous dry-deposition algorithm from GEOS-Chem v10-01 (https://www.geos-chem.org, last access: 4 October 2024) implemented to run in single-point mode with the option to use on-site meteorology and canopy characterizations (i.e., leaf area index (LAI), canopy height, and land type classification). Section S1 in the Supplement details the formulations used in GEOS-Chem for the resistance-in-series components of Eq. (2), which were also detailed in the recent literature (Wong et al., 2019). Briefly, aerodynamic resistance R_a to the turbulent transport of scalars from a reference height z down to the roughness length z_o of the surface is computed following standard surface layer flux–gradient relationships (Wesely and Hicks, 1977). The quasi-laminar boundary layer resistance R_b is estimated following the semiempirical formulation from Wesely and Hicks (1977) and has dependency on both friction velocity u_∗ and species-specific molecular diffusivity D_x. Surface resistance R_c is computed following a big-leaf scheme based on the W89 algorithm, modified for application to the global scale (Wang et al., 1998).

2.1.2 Non-stomatal branch of Z03 dry-deposition algorithm

The Z03 dry-deposition algorithm includes several updates to the W89 scheme. Z03 is used in the Canadian Air and Precipitation Monitoring Network (CAPMoN) (Zhang et al., 2009) as well as in air quality models (e.g., Zhang et al., 2002a) and was recently compared globally to the W89 scheme as implemented in GEOS-Chem for V_d(O₃) (Wong et al., 2019). We implement a stand-alone version of the non-stomatal branch of the Z03 dry-deposition algorithm for NO₂ to enable evaluation against nocturnal eddy-covariance-inferred V_d(NO₂) at Harvard Forest. Leaf stomata are treated as being fully closed at night in the Z03 scheme; accordingly, we treat stomatal resistance as infinite. The Z03 algorithm assigns scale factors (their α and β) for non-stomatal NO₂ uptake relative to inverse surface resistances (conductances) for SO₂ and O₃ of 0 and 0.8, respectively, resulting in parameterized NO₂ deposition velocities 10 %–20 % smaller than for O₃ (Zhang et al., 2002a). The canopy is flagged as wet from dew following the formalism adopted by Brook et al. (1999), with dependence on cloud fraction, temperature, dew point, and u_∗. We estimate the snow cover fraction from snow depth following Zhang et al. (2003a). Following a similar approach to Wu et al. (2018), we compute component surface resistances for the mixed forest as an average of deciduous broadleaf and evergreen needleleaf land-type-specific values from Z03, weighted by the LAI-determined deciduous and coniferous fractions for Harvard Forest of ∼ 80 % and 20 %, respectively (Fig. S5 in the Supplement). As with the stand-alone GEOS-Chem dry-deposition algorithm, we use on-site meteorology and canopy characterizations when available (Sect. 2.2.2).

2.2 Above-canopy dry-deposition velocities inferred from eddy covariance measurements

We evaluate the gaseous dry-deposition scheme from GEOS-Chem against eddy-covariance-inferred deposition velocities over two temperate forests in the USA. First, we compare to deposition velocities from Nguyen et al. (2015) for species found to dry deposit with minimal surface resistance. Being able to neglect the complexities of a surface resistance scheme allows for a more direct evaluation of R_a and R_b components of the resistance-in-series pathway used in the parameterization of V_d (Wu et al., 2021). Second, for an in-depth evaluation of simulated V_d(NO₂) and V_d(NO_y), we use a publicly available long-term hourly dataset of eddy covariance flux observations of NO₂ and NO_y from Harvard Forest, supported with ancillary measurements including NO_y component concentrations, meteorological observations, and canopy characteristics.

2.2.1 Talladega National Forest: H₂O₂, HMHP, and HNO₃

Nguyen et al. (2015) present a novel dataset containing eddy-covariance-inferred deposition velocities of 16 gaseous species, including species found to deposit with negligible surface resistance: H₂O₂, hydroxy methylhydroperoxide (HMHP), and HNO₃. Observations were taken at the Centreville (CTR) Southeastern Aerosol Research and Characterization Study (SEARCH) site (32.90289° N, 87.24968° W) near Brant, Alabama, USA, in June 2013. The CTR site is situated in a grassy clearing in the Talladega National Forest – a mixed forest consisting of coniferous and deciduous tree species, with a mean canopy height h_c of ∼ 10 m and LAI of 4.7 m² m⁻². Eddy covariance flux observations were measured at 22 m a.g.l (meters above ground level). The analysis of Nguyen et al. (2015) includes daytime mean (10:00–15:00 local time (LT)) deposition velocities averaged across 5 ideal days in June 2013 when winds had exclusively forest fetch. To compare with the reported daytime deposition velocities for H₂O₂ (5.2 ± 1.1 cm s⁻¹), HMHP (4.1 ± 1.1 cm s⁻¹), and HNO₃ (3.8 ± 1.3 cm s⁻¹), we average the R_a and R_b components of the stand-alone GEOS-Chem dry-deposition algorithm, applied at the location of the CTR site for the hours of 10:00–15:00 LT on the aforementioned days. The meteorological inputs required to compute the R_a and R_b components of the algorithm ($u_{\ast },T,P$, and sensible heat flux) were obtained from NASA's Goddard Earth Observing System (GEOS) Forward Processed (FP) assimilated meteorological fields (Lucchesi, 2013) at the native horizontal resolution of $\mathrm{1}/\mathrm{4}$° × $\mathrm{5}/\mathrm{16}$°, which Nguyen et al. (2015) note are in excellent agreement with the values measured at the CTR site during this period.

2.2.2 Harvard Forest: NO₂ and NO_y

The utility of the Harvard Forest Environmental Monitoring Site (HFEMS) for evaluating parameterizations of atmosphere–surface exchange stems from the extensive datasets of meteorological and trace gas observations spanning many months to years at high temporal (hourly) resolution. The HFEMS is located in central Massachusetts, USA (42.54° N, 72.18° W; 340 m a.s.l. (meters above sea level)) and is situated in a mature mixed forest (h_c ≃ 20 m), with a summertime LAI and deciduous LAI (DLAI) of 4.3 and 3.4, respectively (Fig. S5). Due to prevailing westerly winds, emissions from Boston (100 km to the east) rarely influence the site. Cool, dry, and unpolluted air from the northwest and warm, moist, anthropogenically influenced air from the southwest are the predominant influences at this site (Horii et al., 2005).

Munger et al. (1996) have described the methodology of the long-term above-canopy (29 m) total nitrogen oxide (NO_y) concentration measurements for eddy covariance flux computation, as well as other details of the HFEMS. Measurements of above-canopy PAN concentrations were added in April 2000 (Horii et al., 2005). Eddy covariance fluxes of NO₂ along with above-canopy (22 m) measurements of HNO₃ concentrations were made at the HFEMS from April through November 2000 (Horii et al., 2004).

Trace gas data from the HFEMS used in this study, specifically hourly eddy covariance fluxes of NO_y and NO₂ and hourly concentrations of NO_y, NO, NO₂, PAN, and HNO₃, are publicly available from the Harvard Forest Data Archive (Munger and Wofsy, 2023). Exchange velocities (V_ex) are computed herein by normalizing reported hourly NO_y and NO₂ eddy covariance fluxes by the respective ambient hourly concentrations. Equating V_ex to V_d assumes that the observed flux is due to surface deposition only. Processes causing deviation from this assumption are discussed in later sections and include surface emission of NO, chemical flux divergence of NO₂, and a potential non-zero canopy accumulation rate of NO_y. Eddy covariance fluxes have fewer errors under conditions where turbulence is well developed (Baldocchi, 2003; Cherin et al., 2015; Goulden et al., 1996; Nguyen et al., 2015). Turbulent threshold u_∗ values in the range of 0.15–0.35 m s⁻¹ (median 0.23 m s⁻¹) have been found to be representative of multiple sites across many years (Cherin et al., 2015). Herein, following the approach of Wu et al. (2011), periods of low surface layer turbulence (u_∗ < 0.2 m s⁻¹) have been omitted from the analysis, resulting in ∼ 25 % of the hourly values of nocturnal V_d(NO₂) and 18 % of the hourly values of V_d(NO_y) being removed from the dataset. Outliers in the remaining hourly V_d(NO₂) and V_d(NO_y) time series were identified via the method of median absolute deviation (MAD) (Leys et al., 2013), where hourly values outside of the median $\pm \mathrm{3}\times$ MAD were removed from calculations of subsequent means; ∼ 20 % of the hourly nocturnal V_d(NO₂) and 10 % of the hourly V_d(NO_y) were removed from the u_∗ filtered dataset. Overall, 60 % of the nocturnal V_d(NO₂) and 74 % of the V_d(NO_y) hourly time series were retained for analysis after application of these turbulence and outlier filters. Figure S6 in the Supplement depicts monthly fractional coverage of hourly measurements of above-canopy trace gas concentrations and eddy-covariance-observed exchange velocities from 2000 to 2002, filtered for conditions of low turbulence.

Meteorological input variables required in the parameterization of V_d were taken from the HFEMS data archive; specifically, P, T, RH, $u_{\ast },$ and sensible heat flux (Munger and Wofsy, 2024) and incoming solar radiation (Fitzjarrald and Sakai, 2023) were available at an hourly temporal resolution throughout the study period. Cloud fraction and snow depth were the only required meteorological variables not available and were instead taken from NASA's Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2)-assimilated meteorological fields (Gelaro et al., 2017). Figure S7 in the Supplement depicts comparisons of hourly observations of u_∗, sensible heat flux, downward shortwave radiation, T, P, and RH made over Harvard Forest to coincident values from MERRA-2-assimilated meteorology, depicting good to excellent agreement. Canopy-specific inputs to the parameterization of V_d include roughness length z_o, displacement height d, LAI, and deposition land type. Values for z_o and d were estimated as $\mathrm{1}/\mathrm{10}$ and $\mathrm{2}/\mathrm{3}$ of canopy height, respectively – values representative of many vegetative surfaces (Garratt, 1992; Oke, 1987), including z_o for Harvard Forest (Wu et al., 2011). We estimate daily LAI values from a spline fit to daily plant area index (PAI) measurements from the HFEMS over April–December for the years 1998–2015 (Matthes et al., 2024), corrected for the reported stem and twig area index (STAI) of 0.9 m² m⁻² noted for this canopy (Horii et al., 2005). Estimated climatological daily LAI values range from ∼ 0.9 m² m⁻² in winter to 4.3 m² m⁻² in summer, in good agreement with MODIS LAI, used by GEOS-Chem, at the location of Harvard Forest (Fig. S5).

2.3 Measured diffusion coefficients of atmospherically relevant molecules

A main result of the work of Nguyen et al. (2015) was the importance of molecular diffusion in atmosphere–surface exchange of rapidly depositing compounds, where it was shown that maximum daytime dry-deposition velocities scale with the inverse square root of molecular mass, as do gas-phase diffusion coefficients (Poling and Prausnitz, 2004). To evaluate the calculation of molecular diffusivities used in the parameterization of gaseous dry-deposition velocities in GEOS-Chem, we conducted a literature search to compile a list of measured diffusion coefficients for atmospherically relevant molecules, consisting of 23 inorganic and 17 organic species (Table S1 in the Supplement). Diffusion coefficients (D) measured in either air or N₂ near standard temperature and pressure (STP) were corrected to STP following Langenberg et al. (2020):

$$\begin{array}{}\text{(3)}& D=D_o\left({\displaystyle \frac{P_o}{P}}\right){\left({\displaystyle \frac{T}{T_o}}\right)}^b,\end{array}$$

where we set the temperature power dependence to b=1.75 following Fuller's method, a semiempirical technique for the estimation of binary gas-phase diffusion coefficients (Fuller et al., 1966), which is discussed further in Sect. 3.2.

2.4 Measurements of surface-specific deposition velocities for NO₂

Surface-specific NO₂ uptake coefficients (${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$) to both foliar and non-foliar forest elements facilitate bottom-up estimates of bulk-canopy R_c(NO₂) and the resulting V_d(NO₂) to forest environments when the corresponding surface area scale factors (i.e., LAI and STAI) and meteorological data are available. From the literature values of surface-specific deposition velocities $V_{\mathrm{d}}^{\mathrm{surf}}$, we infer ${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$ for both non-foliar and foliar materials:

$$\begin{array}{}\text{(4)}& {\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}={\displaystyle \frac{\mathrm{4}{V}_{\mathrm{d}}^{\mathrm{surf}}}{\stackrel{\mathrm{‾}}{v_{\mathrm{t}}}}},\end{array}$$

where $\stackrel{\mathrm{‾}}{v_{\mathrm{t}}}$ is the mean thermal speed of NO₂. Table S2 in the Supplement contains the literature values of $V_{\mathrm{d}}^{\mathrm{surf}}$ and the associated experimental temperatures that were used to infer values of ${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$ to foliar surfaces of deciduous and coniferous species under nocturnal/dark conditions, non-foliar forest materials (bark and forest floor), snow, and fabricated materials. With the exception of deposition to snow, the literature values of $V_{\mathrm{d}}^{\mathrm{surf}}$ are from chamber studies where mechanically mixed chamber air enables the direct estimation of ${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$ through Eq. (4); i.e., turbulent (R_a) and quasi-laminar (R_b) resistances may be neglected for species with slow surface uptake (such as NO₂) since $R_{\mathrm{c}}\gg R_{\mathrm{a}}+R_{\mathrm{b}}$. Values of leaf-level $V_{\mathrm{d}}^{\mathrm{surf}}$ for both deciduous and coniferous species were averaged across periods of minimum stomatal conductance resulting from the absence of photosynthetically active radiation (PAR) and are interpreted herein for the purpose of making the resulting ${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$ value non-stomatal. Table S2 also includes the corresponding surface areas used for flux normalization; care must be taken when comparing surface-specific $V_{\mathrm{d}}^{\mathrm{surf}}$ and ${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$, as various surface areas are used (i.e., planar, geometric, projected leaf area, and total leaf area). Hanson et al. (1989) report $V_{\mathrm{d}}^{\mathrm{surf}}$ to coniferous species normalized to total leaf area, as stomata are distributed across the whole needle surface (amphistomatic). Often, studies normalize to projected leaf area, as is routinely done for deciduous leaves, which generally have stomata on the lower (abaxial) leaf surface. Failing to recognize this difference would result in a misrepresentation of $V_{\mathrm{d}}^{\mathrm{surf}}$ and inferred ${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$ by a factor of ∼ 2.7 for coniferous species (Riederer et al., 1988).

Table 1 summarizes the values of ${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$, including those used in Sect. 3.3.4 to compute nocturnal bottom-up estimates of bulk-canopy R_c(NO₂) and the resulting V_d(NO₂) over Harvard Forest. Also included are the associated surface area scale factors α for which ${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$ is to be applied and the relative humidities over which measurements were made. We suggest that NO₂ uptake to the surfaces listed in Table 1 may result from heterogeneous hydrolysis of NO₂ following Reaction (R1), with variability between surfaces primarily a result of differences in microscopic surface area supporting adsorbed water. Some of the studies measuring foliar uptake of NO₂ under conditions where stomatal aperture should be at a minimum conclude that uptake could occur to the interior of leaves via partially open stomata rather than non-stomatally to the exterior leaf surfaces (Breuninger et al., 2013; Chaparro-Suarez et al., 2011; Delaria et al., 2020; Rondón et al., 1993). Our assumption of nocturnal stomatal closure with deposition of NO₂ to the exterior of leaves via Reaction (R1) is discussed further in Sect. 4. To help contextualize values of ${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$, Sect. S3 in the Supplement provides a brief literature review of uptake coefficients for NO₂ to hydrated surfaces.

Table 1 Surface-specific NO₂ uptake coefficients ${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$ inferred from the literature values of surface-specific deposition velocities^a following Eq. (4). Also included are surface area scale factors α for which ${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$ is to be applied and the surface areas and relative humidities over which surface-specific deposition measurements were made.

^a Surface-specific deposition velocities ($v_{\mathrm{d}}^{\mathrm{surf}}$) were taken from chamber studies, with the exception of uptake to snow, which was measured via the eddy covariance technique. Table S2 in the Supplement contains study-specific details. ^b Surface area used to normalize surface-specific deposition fluxes in the computation of $v_{\mathrm{d}}^{\mathrm{surf}}$ (Table S2). ^c Values used in Sect. 3.3.4 to compute bottom-up estimates of nocturnal bulk-canopy V_d(NO₂) over Harvard Forest. ^d Foliar uptake was measured under conditions of minimal stomatal aperture, i.e., dark conditions. We assume this uptake to be non-stomatal (Sect. 4). ^e Multi-study mean values computed herein (Table S2). ^f Reported $v_{\mathrm{d}}^{\mathrm{surf}}$ normalized to projected leaf areas were scaled herein to reflect uptake to total leaf surface area – a factor of 2.7 for coniferous needles (Riederer et al., 1988). ^g References for surface-specific $V_{\mathrm{d}}^{\mathrm{surf}}({\mathrm{NO}}_{\mathrm{2}}$) are (1) Hanson et al. (1989), (2) Grøntoft and Raychaudhuri (2004), (3) Rondón et al. (1993), (4) Stocker et al. (1995), (5) Wang et al. (2020), (6) Delaria et al. (2020), (7) Delaria et al. (2018), (8) Chaparro-Suarez et al. (2011), and (9) Breuninger et al. (2013). NA denotes not applicable.

3 Measurement–model comparisons and updates

Table 2 summarizes modifications made herein to the offline gaseous dry-deposition parameterization from GEOS-Chem discussed throughout this section. Briefly, parameterization P1 is equivalent to the dry-deposition scheme in GEOS-Chem, which references deposition from grid box centers (GBCs) of the lowest model level (∼ 60 m). Serial modifications to P1 include changes to the height at which dry deposition is referenced (P2), changes to the calculation of aerodynamic resistance (P3) and molecular diffusivities (P4), updating non-stomatal surface resistance for NO₂ following the Z03 scheme (P5) and subsequent replacement with a scheme that represents heterogeneous hydrolysis on deposition surfaces (P6, P7), and finally, implementation of empirical updates to the non-stomatal uptake of PAN (P8). We begin by evaluating parameterizations P1–P4 by comparing to measured dry-deposition velocities from Nguyen et al. (2015), where it was noted that above-canopy deposition velocities for H₂O₂, HMHP, and HNO₃ corresponded to computed theoretical maximums (i.e., $V_{\mathrm{d},\max }\simeq (R_{\mathrm{a}}+R_{\mathrm{b}}{)}^{-\mathrm{1}}$), thus enabling a more direct evaluation of the deposition pathway consisting of resistances R_a and R_b, as discussed in Sect. 3.1 and 3.2, respectively. Parameterizations P5–P7 are evaluated by comparing to both above-canopy nocturnal V_d(NO₂) observed at the HFEMS (Sect. 3.3.3) and bottom-up estimates of nocturnal V_d(NO₂) for Harvard Forest from the literature values of surface-specific deposition velocities $V_{\mathrm{d}}^{\mathrm{surf}}$(NO₂) (Sect. 3.3.4). Parameterization P8 is evaluated in Sect. 3.4 in the context of its effects on simulated V_d(NO_y), including a comparison to the above-canopy diel V_d(NO_y) observed at the HFEMS.

Table 2 Modifications to the offline dry-deposition parameterization tested in this study. Parameterization P1 is equivalent to the gaseous dry-deposition scheme in GEOS-Chem (GC). Modifications to P1 include changes to reference height z_ref (P2), to the formulation of aerodynamic resistance R_a (P3), to molecular diffusivity D (P4), and to non-stomatal surface resistances R_c for NO₂ (P5–P7) and PAN (P8).

^a Dry-deposition reference height is z_GBC ≃ 60 m, z_TNF=22 m, and z_HFEMS=29 m. ^b Mean free path (mfp) is held constant across depositing gases. ^c Empirical fit of non-stomatal cuticular deposition (Turnipseed et al., 2006) modified herein for LAI (Sect. S5 in the Supplement).

3.1 Updates to the calculation of aerodynamic resistance

Table 3 contains an evaluation of simulated atmospheric resistances (R_a and R_b) in parameterizations P1–P4 by comparing to measured daytime deposition velocities for rapidly depositing species from Nguyen et al. (2015). In contrast to the findings of Wu et al. (2021) that show excellent model–measurement agreement of peak daytime V_d(HNO₃) and V_d(H₂O₂) between the Z03 scheme and the dataset from Nguyen et al. (2015), parameterization P1 overestimates daytime mean deposition velocities for the following rapidly depositing species: H₂O₂ (+15 %), HMHP (+41 %), and HNO₃ (+52 %). Nguyen et al. (2015) found excellent agreement between hourly GEOS-FP-assimilated meteorology at this site (used herein for computation of V_d in Table 3) and measured values, including u_∗ and sensible and latent heat fluxes. GEOS-FP fields report a summertime z_o=2.2 m for the 0.25° × 0.3125° grid cell that includes the CTR site – greater than would be expected at this site given the local 10 m canopy height. Prescribing z_o to be 10 % of h_c in parameterization P1b, in accordance with conventionally used values for natural vegetation and in agreement with an updated land use module developed for GEOS-Chem (Geddes et al., 2016), results in a 35 % increase in R_a and notable reductions in V_d high biases. However, following the computation of R_a in GEOS-Chem, P1b computes R_a from a reference height z of ≃ 60 m despite a measurement height of 22 m at the CTR site, while neglecting a displacement height d. Neglecting d from the computation of R_a in Eq. (S1) in the Supplement increases daytime R_a in parameterization P1b by 1 % when referenced from 60 m and 9 % when referenced from 22 m. Although the greatest sensitivity of R_a to z occurs in proximity to z_o (Fig. S2), the difference between R_a computed from an above-canopy measurement height and the typical heights from which global CTMs reference dry deposition can be significant (Figs. S1 and S2). Referencing R_a from the CTR measurement height of 22 m in parameterization P2 results in a 23 % decrease in daytime R_a, returning V_d biases to P1 levels (Table 3).

Table 3 Effects of updates to the calculation of aerodynamic resistance (R_a) and quasi-laminar sublayer resistance (R_b) on simulated daytime (10:00–15:00 LT) dry-deposition velocities (V_d) over Talladega National Forest (temperate, mixed) for three rapidly depositing species. Serial modifications to the base parameterization P1 are highlighted, i.e., PX (update). Shown are the mean quantities ± standard deviations about the hourly time series^a and the normalized mean bias (NMB) between simulated and measured^b V_d.

^a Mean quantities are averaged across the five daytime periods in June 2013 that Nguyen et al. (2015) used in their analysis of eddy-covariance-observed deposition velocities (Sect. 2.2.1). ^b Measured (eddy covariance) daytime (10:00–15:00 LT) V_d(H₂O₂)=5.2 ± 1.1 cm s⁻¹, V_d(HMHP) = 4.1 ± 1.1 cm s⁻¹, and V_d(HNO₃)=3.8 ± 1.3 cm s⁻¹ (Nguyen et al., 2015). ^c Table 2 contains a list of parameterization updates. Surface resistance (R_c) was set to 1 s m⁻¹ following the minimum allowed in GEOS-Chem (Sect. S1 in the Supplement). ^d Roughness length (z_o) was set to 10 % of canopy height (h_c) for parameterizations P1b–P8.

Considering the fact that the CTR and HFEMS measurement heights, ∼ 2 h_c and 1.5 h_c, respectively, are at the upper limits of the roughness sublayer (RSL), a region where turbulent mixing in the wake of roughness elements is enhanced by a factor of 2 to 3 above that predicted by Monin–Obukhov (M–O) similarity theory (see Sect. S2.4 in the Supplement for a review of RSL mixing), R_a computed according to M–O similarity theory following Eq. (S1) may be a slight underestimate due to non-zero horizontal winds at z_o resulting from enhanced downward mixing of momentum. To quantify this effect, parameterization P3 computes R_a corrected for RSL mixing, which allows u(z_o) > 0 m s⁻¹ (Eq. S13), resulting in a small (5 %) increase in R_a at the CTR measurement height under the daytime conditions of Table 3 and even smaller changes to V_d given the influence of R_b (Sect. 3.2). As demonstrated herein and in agreement with previous work (Simpson et al., 1998), it may be appropriate to neglect the effects of RSL on depositing species when referenced from a height of at least 1.5–2 h_c. However, studies endeavoring to understand bidirectional exchange or the dispersion of near-surface emissions should consider the effect of asymmetrical R_a that the RSL imposes (Sect. S2.2 and S2.4).

3.2 Updates to the calculation of molecular diffusivities

As seen in Table 3, updates to the calculation of R_a failed to address high biases in simulated deposition velocities of rapidly depositing species. The higher-molecular-weight species HMHP and HNO₃ exhibit a greater high bias in V_d, 36 % and 47 %, respectively, than the lower-molecular-weight species H₂O₂ (11 %). Given the dependence of V_d,max on molecular diffusivity through its influence on R_b (Eq. S2) (Meyers et al., 1989), we evaluate the calculation of molecular diffusion coefficients in GEOS-Chem against measured values for atmospherically relevant molecules.

Figure 1 depicts a large high bias in calculated diffusion coefficients from the dry-deposition module of GEOS-Chem, which uses the Chapman–Enskog theory for binary diffusivity (Seinfeld, 1986). The bias results from the use of a constant collision diameter of σ=2.7 Å with air for all species – an underestimate for many atmospherically relevant molecules, i.e., σ for O₃ with air is 3.793 Å (Massman, 1998; Poling and Prausnitz, 2004). The collision diameter σ is a pairwise characteristic length scale of the Lennard–Jones intermolecular force, which is not readily available for many atmospheric trace gases (Tang et al., 2014). Several semiempirical methods have been proposed for the estimation of D (the diffusion coefficient) in low-pressure binary systems (Poling and Prausnitz, 2004). Fuller et al. (1966) developed a simple and generalized semiempirical correlation equation for the estimation of binary gas-phase diffusion coefficients using additive atomic diffusion volumes V_i for each species $\sum _A{V}_i$ and $\sum _B{V}_i$. The diffusion coefficient D [cm² s⁻¹] for trace gas A in bath gas B is given by

$$\begin{array}{}\text{(5)}& D={\displaystyle \frac{{\mathrm{10}}^{-\mathrm{3}}{T}^{\mathrm{1.75}}{\left(\mathrm{1}/M_A+\mathrm{1}/M_B\right)}^{\mathrm{1}/\mathrm{2}}}{P{\left[{\left(\sum _A{V}_i\right)}^{\mathrm{1}/\mathrm{3}}+{\left(\sum _B{V}_i\right)}^{\mathrm{1}/\mathrm{3}}\right]}^{\mathrm{2}}}},\end{array}$$

where P is the pressure [atm], T is the temperature [K], and M is the molecular mass [g mol⁻¹]. Atomic and, in some cases, molecular diffusion volumes were obtained from regression analysis of 153 binary systems across 340 T–P states and are summarized in Poling and Prausnitz (2004), Tang et al. (2014), and Tang et al. (2015). As seen in Fig. 1, diffusion coefficients computed using Fuller's method result in a much-improved comparison to measurements, with better agreement to organic species (R²=0.99 and NMB = −3 %) than to inorganics (R²=0.88 and NMB = 13 %), consistent with the findings of Tang et al. (2014, 2015).

Figure 1 Measured diffusion coefficients of atmospherically relevant molecules in air or N₂ at STP are compared to the calculated values. Molecular diffusivities calculated following the method used in GEOS-Chem (P1–P3) are compared to those calculated following Fuller's method and Graham's law (referenced from $D_{{\mathrm{CO}}_{\mathrm{2}}}$). Measured and computed (Fuller's method) values are listed in Table S1 in the Supplement.

Figure 1 also depicts molecular diffusion coefficients approximated by Graham's law of effusion, i.e., $D_{\mathrm{1}k}=D_{\mathrm{2}k}\sqrt{M_{\mathrm{2}}/M_{\mathrm{1}}}$ (Mason and Evans, 1969), where (continuum) diffusion coefficients are approximated by the Knudsen diffusion coefficients D_k – an oversimplification of Eq. (5) and a strategy commonly used in the atmospheric science community (Nguyen et al., 2015; Weber and Renenberg, 1996; Wesely, 1989). The resulting diffusion coefficients scaled from measured $D_{{\mathrm{CO}}_{\mathrm{2}}}$ correlate well with the measured values (R²=0.91), with normalized mean biases (NMBs) to inorganic and organic species of 3 % and 20 %, respectively. In a review of molecular diffusivities of atmospherically relevant molecules, Massman (1998) note that misapplication of Graham's law to molecular diffusivities can lead to errors of up to 23 %. Measured and computed diffusion coefficients from Fuller's method, assuming air as the bath gas, are presented in Fig. 1 and tabulated in Table S1. We do not differentiate between diffusivity measurements carried out in air and in N₂, as differences are expected to be small, i.e., a 2 % difference in $D_{{\mathrm{O}}_{\mathrm{3}}}$ at STP in air vs. in N₂.

Parameterization P4 computes R_b using measured diffusion coefficients when available and diffusion coefficients according to Fuller's method in the absence of measured values. Diffusion coefficients are adjusted to ambient T-P following Eq. (3) prior to calculating R_b. Eliminating the high bias in calculated molecular diffusivities resulted in a near doubling of R_b for the species in Table 3 and a much-improved comparison to the daytime deposition velocities for the higher-molecular-weight species HMHP (NMB −1 %) and HNO₃ (NMB 4 %). The increase in R_b for H₂O₂ results in a low bias of −15 %, which is well within the large relative uncertainty for R_b due to variations in canopy structures (Massman et al., 1994; Sievering et al., 2001).

Molecular diffusivity is also involved in the calculation of R_c via its influence on stomatal resistance r_s, which is scaled by the ratio $D_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}/D_x$ in the dry-deposition parameterizations commonly used in chemical transport models (Wesely, 1989; Zhang et al., 2003a). The effect of updated molecular diffusivity on R_c in GEOS-Chem is significant for molecules that dry deposit under stomatal control, i.e., species with low aqueous solubility or surface reactivity, and is discussed in Sect. S5 in the Supplement.

3.3 Nocturnal dry deposition of NO₂ over Harvard Forest

3.3.1 Eddy-covariance-inferred V_d(NO₂)

Nocturnal hourly eddy covariance NO₂ fluxes and the resulting exchange velocities V_ex(NO₂) over Harvard Forest from April to November 2000 are shown in Fig. 2 as a function of NO₂ concentration. We restrict our analysis to nighttime (20:00–04:00 local standard time (LST)), when above-canopy NO₂ : NO_x is ∼ 1, and photochemical flux divergence of the NO–NO₂–O₃ triad due to the presence of a vertical gradient in irradiance through the forest canopy (Gao et al., 1993) is absent. As seen in the top panel of Fig. 2, nocturnal fluxes of NO₂ over Harvard Forest are predominantly (∼ 70 %) downward, especially at higher ambient NO₂ concentrations. Nocturnal mean (median) ±1σ fluxes of NO₂ from April to November are −0.8 (−0.3) ± 2 ppb cm s⁻¹. These downward (p < 0.01) above-canopy aggregate fluxes of NO₂ are comparable in magnitude to counteracting summertime nocturnal soil NO emissions, estimated by Munger et al. (1996) through a mass–balance approach to be ∼ 0.9 µmol m⁻² h⁻¹ (3.5 ng N m⁻² s⁻¹ or 0.62 ppb cm s⁻¹) at the HFEMS. Munger et al. (1996) note that nocturnal NO is elevated near the forest floor, and Horii et al. (2004) find decreasing within-canopy nocturnal NO profiles at Harvard Forest, with above-canopy concentrations and fluxes indistinguishable from zero despite net downward fluxes of NO₂, presumably due to the titration of soil-emitted NO by O₃ on a timescale much shorter (minutes) than in-canopy vertical mixing, followed by nocturnal canopy loss processes for NO₂. Previous studies have noted the importance of knowledge of local soil NO emissions and within-canopy processes involving NO_x when interpreting above-canopy NO₂ fluxes (Delaria and Cohen, 2020; Eugster and Hesterberg, 1996; Flechard et al., 2011; Min et al., 2014).

In an effort to isolate the contribution that dry deposition makes to the above-canopy nocturnal eddy covariance fluxes (F_EC) of NO₂, we infer V_d(NO₂) following Eq. (6) to account for the effects of nocturnal chemical flux divergence (V_chem) and to counteract the soil NO emissions that are assumed to rapidly titrate with O₃ and ventilate the canopy as NO₂ (F_soil). The resulting observation-inferred V_d(NO₂) is a best estimate of the nocturnal dry-deposition pathway with which to evaluate parameterizations:

$$\begin{array}{}\text{(6)}& V_{\mathrm{d}}+V_{\mathrm{chem}}=-V_{\mathrm{ex}}=-{\displaystyle \frac{\left(F_{\mathrm{EC}}-F_{\mathrm{soil}}\right)}{\left[{\mathrm{NO}}_{\mathrm{2}}\right]}},\end{array}$$

where V_ex is the eddy-covariance-observed NO₂ exchange velocity, which does not assume predominant deposition and therefore has sign convention analogous to F_EC; V_chem represents an estimate of below-sensor nocturnal chemical loss of NO₂ via formation and loss of N₂O₅, limited by the rate of oxidation of NO₂ with O₃ (Browne and Cohen, 2012; Jacob, 2000). We use an estimate of the maximum rate of nocturnal chemical loss of NO₂ proposed by Horii (2002) in their analysis of the dataset used herein, V_chem ≃ 0.05 cm s⁻¹, which translates to a below-sensor (< 29 m) nocturnal chemical lifetime of NO₂ to oxidation by O₃ of ∼ 16 h. The bottom panel of Fig. 2 includes hourly values of V_ex(NO₂), both uncorrected and corrected for soil NO. Values of F_soil used in Eq. (6) are less than the peak summertime forest floor estimate from Munger et al. (1996) due to seasonality and within-canopy loss processes. Hourly estimates of F_soil were calculated by scaling the reported summertime nocturnal soil NO emission flux at Harvard Forest, F_NO, summer (0.62 ppb cm s⁻¹), by GEOS-Chem-simulated seasonality κ and a parameterized canopy reduction factor (CRF):

$$\begin{array}{}\text{(7)}& F_{\mathrm{soil}}\left(\mathrm{h}\right)=F_{\mathrm{NO},\mathrm{summer}}\mathit{\kappa }\left(\text{month}\right)\left[\mathrm{1}-\text{CRF}\left(\mathrm{h}\right)\right].\end{array}$$

Month-specific κ scale factors were obtained by normalizing simulated monthly mean nocturnal soil NO emission, which was output at the location of the HFEMS from a high-resolution (0.25° × 0.3125°) GEOS-Chem simulation, by the peak monthly mean simulated emission (July at the location of HFEMS). GEOS-Chem simulated soil NO emission in the region of Harvard Forest exhibits significant seasonality, with a winter minimum that is a small fraction (< 5 %) of the summertime maximum (Fig. S3). Section S4 in the Supplement describes the parameterization of CRF used in GEOS-Chem and in Eq. (7). As NO₂ surface uptake resistance R_c(NO₂) is used in the calculation of CRF, values are parameterization specific (Fig. S3), with larger values of CRF resulting from lower values of R_c(NO₂).

Figure 2 (a) Nocturnal (20:00–04:00 LST) hourly eddy covariance NO₂ fluxes and (b) the resulting exchange velocities V_ex(NO₂) as a function of NO₂ concentration. These publicly available measurements (Sect. 2.2.2) were taken over an established mixed forest (Harvard Forest) from April to November 2000. Estimated soil NO flux (assumed to ventilate the canopy as NO₂) was subtracted from measured hourly NO₂ fluxes in order to estimate V_ex(NO₂) corrected for soil NO (Eqs. 6–7). Included in (a) are linear fits and the associated 95 % prediction intervals. Dashed lines in (b) depict boundaries of an outlier filter applied to hourly V_ex(NO₂) prior to the calculation of means (Sect. 2.2.2). Data points excluded from the analysis based on visual inspection are circled. Hourly observations made under conditions of low turbulence (u_∗ < 0.2 m s⁻¹) were excluded from this analysis.

Figure 3 Observation-inferred (eddy covariance) and simulated monthly mean nocturnal (20:00–04:00 LST) NO₂ deposition velocities V_d(NO₂) over Harvard Forest. Simulated values are coincidently sampled with hourly observations prior to averaging. Also depicted is observation-inferred V_d(NO₂) corrected for soil NO emission using simulated soil NO canopy reduction factors (CRFs) from parameterizations P4, P6, and P7. Standard deviations about simulated monthly mean values, as well as observation-inferred monthly mean values uncorrected for soil NO, are depicted as shaded areas. Month-specific ranges in number of hourly observations used in the calculation of monthly means are indicated as vertical lines and result from parameterization-specific soil NO corrections causing differential outlier filter exceedance (Fig. 2 and Table 4). Insufficient data prevented analysis for May and June (Fig. S6).

Figure 3 depicts monthly nocturnal V_d(NO₂) inferred from Eq. (6) over the HFEMS from April to November 2000 alongside coincidently sampled simulated values from parameterizations P4 to P7. Table 4 presents observation-inferred and simulated values of V_d(NO₂) aggregated across all months, as well as associated NO₂ lifetimes to dry deposition from the 29 m measurement height. We begin the discussion of eddy-covariance-inferred bulk-canopy V_d(NO₂) below, followed by discussions in Sect. 3.3.2 of (i) the large low bias in simulated values stemming from the widely used W89 parameterization of surface resistances and (ii) the reduced bias using the Z03 scheme. In Sect. 3.3.3, we evaluate a simple representation of non-stomatal NO₂ uptake following Reaction (R1) against eddy-covariance-inferred V_d(NO₂).

Table 4 Nocturnal (20:00–04:00 LST) NO₂ deposition velocities over Harvard Forest aggregated from April to November. Observation-inferred values with and without estimated soil NO corrections using simulated canopy reduction factors (CRFs) corresponding to parameterizations P4, P6, and P7 are shown along with coincidently sampled simulated values. Measurements under conditions of low turbulence (friction velocity u_∗ < 0.2 m s⁻¹) were excluded from the analysis, leaving 230 hourly observations in the time series (Sect. 2.2.2).

^a Observation-inferred hourly values of V_ex(NO₂), both uncorrected and corrected for soil NO, were subjected to an outlier filter (Fig. 2) prior to the calculation of the arithmetic mean V_d(NO₂); the remaining number of hourly observations are included in brackets adjacent to the corresponding mean values and standard deviations. ^b The ratio of mean NO₂ flux to mean NO₂ concentration ($\stackrel{\mathrm{‾}}{F}/\stackrel{\mathrm{‾}}{\left[{\mathrm{NO}}_{\mathrm{2}}\right]}$) [cm s⁻¹] is included for comparison (Sect. 3.3.1). ^c Calculated at 29 m (measurement height) from aggregate mean deposition velocities. ^d See Table 2 for serial updates. Briefly, P4 has updates to aerodynamic resistance (R_a) and molecular diffusivity (D). P5 computes the resistance to non-stomatal surface uptake of NO₂ following Zhang et al. (2003a). P6 computes the resistance to non-stomatal surface uptake of NO₂ following Reaction (R1), with surface area scale factor α=1 in Eq. (8). P7 is analogous to P6 but with α=2. N/A denotes not applicable.

As previously mentioned, hourly values of observed V_ex(NO₂) were subjected to an outlier filter (Fig. 2) prior to computation of mean values, whereas median and “mean flux to mean concentration” ratios ($\stackrel{\mathrm{‾}}{F}/\stackrel{\mathrm{‾}}{\left[{\mathrm{NO}}_{\mathrm{2}}\right]}$) included in Table 4 were not and instead were computed directly from u_∗-filtered hourly data, as the latter two statistics are less influenced by outliers than arithmetic means. Aggregate values of $\stackrel{\mathrm{‾}}{F}/\stackrel{\mathrm{‾}}{\left[{\mathrm{NO}}_{\mathrm{2}}\right]}$ in Table 4 are in the same units as V_d(NO₂) [cm s⁻¹] and include corrections for F_soil and V_chem, as do the mean and median quantities computed from hourly values following Eq. (6). Assuming a first-order dependence of NO₂ dry deposition on concentration (Eq. 1), computing values of $\stackrel{\mathrm{‾}}{F}/\stackrel{\mathrm{‾}}{\left[{\mathrm{NO}}_{\mathrm{2}}\right]}$ over long averaging times is a strategy to reduce the influence of random variability in deposition velocity estimates, especially under low-NO₂ conditions, as are evident in Fig. 2.

As seen in Fig. 3, monthly mean values of observation-inferred V_d(NO₂) uncorrected for the influence of soil-emitted NO are in the range of 0.1–0.2 cm s⁻¹. Although the variability is large, with standard deviations greater than the mean values, corrections for soil NO (venting from the canopy as NO₂) result in a significant (p < 0.03) increase in nocturnal V_d(NO₂), yielding monthly mean values in the approximate range of 0.2–0.3 cm s⁻¹. Both uncorrected and soil-NO-corrected nocturnal V_d(NO₂) lack discernible seasonality. In Sect. 3.3.4, bottom-up estimates of nocturnal V_d(NO₂) for Harvard Forest are developed in an effort to understand the apparent lack of seasonality in top-down observations. Large variability in the eddy-covariance-observed NO₂ flux and the resulting deposition velocities has been noted in other studies (Eugster and Hesterberg, 1996; Farmer et al., 2006; Geddes and Murphy, 2014), wherein authors restrict their analysis to average values in order to reduce the variability in these complex ecosystem-scale observations (Baldocchi, 2003). Herein, we restrict our analysis to the average values over at least 1 month.

3.3.2 Evaluation of nocturnal V_d(NO₂) from GEOS-Chem and Z03

Parameterization P4, which computes NO₂ surface uptake resistance R_c(NO₂) following the W89 representation in GEOS-Chem, yields a simulated nocturnal V_d(NO₂) that has a nearly 4-fold low bias compared to observations uncorrected for soil NO, increasing to a 7-fold low bias after correcting for soil NO from parameterization P4 (Table 4). This underestimate is driven by the large nocturnal R_c(NO₂) of ∼ 2700 s m⁻¹ in parameterization P4 (Fig. S8), which has been noted in previous studies comparing the W89 algorithm to eddy covariance observations over forest (Horii et al., 2004) and grassland (Eugster and Hesterberg, 1996) ecosystems. Wesely et al. (1982) reported a nocturnal eddy-covariance-observed V_d(NO₂) of 0.05 cm s⁻¹ over a summertime soybean field, similar to the P4 value in Table 4. The authors acknowledge that counteracting soil NO emissions may have resulted in low measured values of above-canopy NO₂ deposition. In their analysis of eddy covariance fluxes of NO₂ over a managed grassland in central Switzerland, Eugster and Hesterberg (1996) found that accounting for counteracting fluxes of soil-emitted NO, oxidized to NO₂ below the height of the sensor (∼ 2.7 to 3.6 ng N m⁻² s⁻¹), resulted in an increase in inferred nocturnal V_d(NO₂) by up to a factor of 2, corresponding to an inferred median value for nocturnal non-stomatal R_c(NO₂) of 700 s m⁻¹ (range 500–950 s m⁻¹) – a surface resistance on the order of 4 times lower than was predicted by the W89 algorithm.

In parameterization P5, R_c(NO₂) is computed following the Z03 scheme, resulting in large increases in simulated nocturnal V_d(NO₂) across all months, as well as increased intra- and inter-month variability (Fig. 3). The Z03 scheme increases non-stomatal NO₂ uptake relative to in the W89 scheme. The representation in the W89 scheme results in stomatal control over uptake, with very low non-stomatal uptake (Sect. S5, Fig. S4). The diurnal behavior of observed V_d(NO₂) relative to V_d(O₃) reported in Wesely et al. (1982) provided support for the use of W89. For use in algorithmic developments to gaseous dry deposition, Zhang et al. (2002a) compiled an updated list of half-reaction redox potentials for species of interest in dry-deposition models, noting greater oxidizing capacity for both NO₂ and NO over their assignments in W89. Zhang et al. (2002a) justify setting the non-stomatal conductance of NO₂ at 80 % that of O₃ – greater than the 10 % used in W89 – based on noted oxidizing capacity and contemporary field studies (Eugster and Hesterberg, 1996; Pilegaard et al., 1998; Rondón et al., 1993; Walton et al., 1997). As seen in Fig. 3, we find the Z03 scheme in parameterization P5 to be in good agreement with eddy-covariance-inferred nocturnal V_d(NO₂) at Harvard Forest, with monthly mean values between the observation-inferred values that are uncorrected and corrected for emissions of soil NO. As seen in Table 4, parameterization P5 compares well to observation-inferred V_d(NO₂) corrected for soil NO with CRF at greater levels of uptake (i.e., P7 levels).

Although nocturnal V_d(NO₂) from parameterization P5 is in good agreement with observation-inferred values over Harvard Forest, it is difficult to justify the increase in parameterized non-stomatal uptake of NO₂ in Zhang et al. (2002a), a predecessor to the Z03 scheme, from half-reaction redox potentials alone. Support for the level of non-stomatal NO₂ uptake in the Z03 scheme comes largely from field observations – many of which are referenced herein in subsequent sections, as we discuss the plausibility of non-stomatal NO₂ uptake being a result of heterogeneous hydrolysis – rather than from suspected NO₂ reduction reactions based on half-reaction redox potentials. Exploring this potential misrepresentation of non-stomatal NO₂ uptake may have important implications for future representation of HONO surface sources in atmospheric CTMs.

3.3.3 Updates to parameterized V_d(NO₂) by representing NO₂ hydrolysis on deposition surfaces

Horii et al. (2004) note that observed nocturnal dry deposition of NO₂ may result from a surface hydrolysis reaction following Reaction (R1). In parameterizations P6–P8, we replace the non-stomatal components of the bulk-surface resistance scheme for NO₂ with a dry-deposition pathway representing NO₂ hydrolysis r_hyd [s m⁻¹], formulated as a collision-limited heterogeneous reaction with ground surfaces (Cano-Ruiz et al., 1993):

$$\begin{array}{}\text{(8)}& r_{\mathrm{hyd}}={\displaystyle \frac{\mathrm{4}}{{\mathit{\gamma }}_{\mathrm{g},{\mathrm{NO}}_{\mathrm{2}}}\stackrel{\mathrm{‾}}{v_{\mathrm{t}}}\mathit{\alpha }}},\end{array}$$

where ${\mathit{\gamma }}_{\mathrm{g},{\mathrm{NO}}_{\mathrm{2}}}$ is a ground uptake coefficient for NO₂ resulting from heterogeneous hydrolysis on deposition surfaces, $\stackrel{\mathrm{‾}}{v_{\mathrm{t}}}$ the mean thermal speed of NO₂, and α a dimensionless scale factor introduced herein to facilitate application of Eq. (8) across land types of varying surface area densities. Lammel and Cape (1996) recommend that Reaction (R1) be parameterized in atmospheric chemistry models using field-derived uptake coefficients, as realistic conditions are difficult to reproduce in the lab. We use the field-derived ground uptake coefficient ${\mathit{\gamma }}_{\mathrm{g},{\mathrm{NO}}_{\mathrm{2}}}$ for Reaction (R1) from VandenBoer et al. (2013), determined from the production of HONO in a winter nocturnal boundary layer in an agricultural region of Colorado, USA. Consistent with the heterogeneous hydrolysis of NO₂ requiring adsorbed water to proceed, VandenBoer et al. (2013) parameterized ${\mathit{\gamma }}_{\mathrm{g},{\mathrm{NO}}_{\mathrm{2}}}$ as a function of RH [ %] according to Eq. (9) to capture the factor of 2 variability in ${\mathit{\gamma }}_{\mathrm{g},{\mathrm{NO}}_{\mathrm{2}}}$ on either side of their best-fit value ($\mathrm{8}\times {\mathrm{10}}^{-\mathrm{6}}$):

$$\begin{array}{}\text{(9)}& {\mathit{\gamma }}_{\mathrm{g},{\mathrm{NO}}_{\mathrm{2}}}={\displaystyle \frac{\mathrm{RH}}{\mathrm{50}}}\mathrm{8}\times {\mathrm{10}}^{-\mathrm{6}}.\end{array}$$

Parameterization P6 computes r_hyd with α=1, resulting in a simulated nocturnal aggregate mean V_d(NO₂) of 0.12 ± 0.02 cm s⁻¹ – a 3-fold increase over P4 and satisfactory agreement with observation-inferred V_d(NO₂) uncorrected for soil-emitted NO (Table 4). However, this is an underestimate by ∼ 50 % when soil NO emissions are accounted for. The larger nocturnal mean CRF of 47 % for parameterization P6 is due to reduced nocturnal R_c(NO₂) (median value ∼ 750 s m⁻¹, Fig. S8), resulting in a small (11 %) decrease in observation-inferred V_d(NO₂) corrected for soil NO (Table 4). Increasing the rate of non-stomatal uptake of NO₂ by computing r_hyd with α=2 in parameterization P7 resulted in a simulated nocturnal aggregate mean V_d(NO₂) of 0.24 ± 0.04 cm s⁻¹, which compares well to the observation-inferred V_d(NO₂) of 0.23 ± 0.35 cm s⁻¹ after correction for soil-emitted NO using a 59 % CRF from P7 (Table 4). Although parameterization P7 results in satisfactory simulation of nocturnal V_d(NO₂) at the HFEMS when averaged across all months, intra- and inter-month variability in observation-inferred V_d(NO₂) is not captured in simulated values (Fig. 3).

Physical justification for the scale factor value α > 1 being necessary to reduce the bias between simulated nocturnal V_d(NO₂) and observation-inferred values corrected for estimated soil NO could stem from the larger surface area available for NO₂ heterogeneous hydrolysis in a mature forest environment compared to the US Midwest wintertime agricultural region over which VandenBoer et al. (2013) derived ${\mathit{\gamma }}_{\mathrm{g},{\mathrm{NO}}_{\mathrm{2}}}$. Heterogeneous reactions not limited by transport or diffusion to reaction surfaces are governed by a collision-limited rate that scales linearly with the surface area to volume ratio of the reaction vessel or environment (Jacob, 2000). Heterogeneous hydrolysis of NO₂ may proceed on any surface accommodating adsorbed water, including foliar surfaces, bark, or elements of the forest floor (i.e., rock, soil, and debris). Despite the hydrophobic nature of many foliar surfaces, thin aqueous films have been observed on coniferous needles (Altimir et al., 2006; Burkhardt and Eiden, 1994) and stomata-bearing surfaces of deciduous leaves (Burkhardt et al., 1999) at ambient humidities well-below saturation. In addition to radiative cooling, elevated humidity within the thin laminar boundary layer surrounding leaves may result from stomatal transpiration (Burkhardt and Hunsche, 2013) and to a lesser extent the hydraulic activation of stomata (HAS) (Burkhardt, 2010), a process discussed further in Sect. 4. Surface area indices [m² m⁻²] for forest components at the HFEMS have been estimated (Fig. S5), including for stems and twigs (STAI = 0.9), coniferous needles (CAI = 0.8), and deciduous leaves (DLAI = 3.4 summertime maximum). Assuming round stems and twigs (Sörgel et al., 2011) and oblate coniferous needles (Oren et al., 1986; Riederer et al., 1988), the total wintertime canopy surface area is estimated as πSTAI + 2.7CAI ≃ 5 m² m⁻². We estimate the summertime canopy surface area to be ∼ 12 m² m⁻² accounting for both sides of deciduous leaves or ∼ 9 m² m⁻² neglecting the non-stomatal adaxial (top) surface of deciduous leaves, in agreement with typical macroscopic surface area indices for temperate and boreal forest canopies (12 m² m⁻², range of 5–14 m² m⁻²) (Lammel, 1999). The surface area of the forest floor, including debris, would also be much greater than the planar ground area and that of tree bark is greater than the simple geometric surface area (Sect. 3.3.4).

The lack of seasonality in observation-inferred nocturnal V_d(NO₂) depicted in Fig. 3 may reflect an inter-seasonal buffering of available surface area for the reaction of above-canopy NO₂ due to increased air parcel mixing throughout the lower canopy in the absence of deciduous leaves (see Sect. 3.3.4). We did not attempt to parameterize non-stomatal deposition of NO₂ to upper- and lower-canopy elements separately in our top-down sensitivity analysis of R_c(NO₂), as is currently the approach in the W89 and Z03 dry-deposition schemes. Due to the lack of discernible seasonal variability in observation-inferred nocturnal V_d(NO₂), above-canopy observations were insufficient to justify the additional variables. We acknowledge that the nocturnal canopy environment is under reduced turbulent mixing compared to daytime conditions when the forest experiences enhanced vertical exchange (Bannister et al., 2023; Sörgel et al., 2011; Thomas and Foken, 2007). Although the daytime surface area available to above-canopy deposition is therefore likely to be greater than that at night, the nighttime sensitivity of V_d(NO₂) to α is much greater than during the day when stomata are open and foliar uptake of NO₂ is a more substantial pathway to deposition than non-stomatal uptake (Fig. S8). Increasing α from 1 to 2 results in a 100 % increase in simulated V_d(NO₂) at night (Table 4) but only a 10 % increase during midday (Sect. S5).

The canopy compensation point for NO₂ is the ambient above-canopy concentration at which point consumption (i.e., dry deposition) and production (i.e., soil emission) are in balance (Duyzer et al., 1995). Studies of above-canopy NO₂ exchange have observed aggregate fluxes to be upward (Min et al., 2014; Vaughan et al., 2016), downward (Coe and Gallagher, 1992; Horii et al., 2004; Walton et al., 1997), and not significantly different from zero (Geddes and Murphy, 2014) – highlighting the importance of knowledge of below-canopy NO_x emission and subsequent uptake and reaction prior to attempting interpretation of above-canopy fluxes. Although foliar compensation points for NO₂ – the concentration below which vegetation is proposed to become a net source of NO₂ – have been observed in the past via leaf-level chamber measurements to be generally < 2 ppb (Geßler et al., 2002; Sparks et al., 2001; Weber and Renenberg, 1996), recent chamber studies using highly specific NO₂ detection methods have failed to observe foliar emission (Breuninger et al., 2013; Chaparro-Suarez et al., 2011; Delaria et al., 2018, 2020; Wang et al., 2020). By restricting analysis to (i) nocturnal conditions when RH is generally high (Fig. S7) and stomata are assumed to be closed (Sect. 4), and to (ii) well-established turbulence (u_∗ > 0.2 m s⁻¹), we find that monthly aggregate V_d(NO₂) is relatively constant in the range of 0.2–0.3 cm s⁻¹ across April–November, with expected large variability on finer timescales (i.e., σ=0.35 cm s⁻¹ across the hourly dataset). Linear regression of hourly nocturnal NO₂ flux vs. ambient NO₂ concentration (Fig. 2) yields a V_ex(NO₂) of −0.40 cm s⁻¹ (p < 0.01) over the entire NO₂ concentration range (up to ∼ 30 ppb) and −0.26 cm s⁻¹ (p < 0.01) when the four outlying hourly observations beyond 20 ppb NO₂ are excluded – consistent with the findings of Horii et al. (2004). An inferred V_d(NO₂) of 0.21 cm s⁻¹ is obtained after subtraction of V_chem=0.05 cm s⁻¹ – similar to the aggregate values presented in Table 4. The y-axis intercept of 0.11 ppb cm s⁻¹ in Fig. 2, although not significant (p > 0.1), is in line with the estimated mean (April–November) above-canopy NO₂ flux of 0.13 ppb cm s⁻¹ resulting from soil NO emission and an average CRF of 59 % from parameterization P7. An empirical CRF of ∼ 70 % is obtained from the ratio of the y-axis intercept (0.11 ppb cm s⁻¹, Fig. 2) to the seasonal mean below-canopy soil NO flux (0.39 ppb cm s⁻¹). An NO₂ canopy compensation point for Harvard Forest, likely due to soil NO emission, is approximated by the x-axis intercept of ∼ 0.4 ppb in Fig. 2.

Given the dependence of r_hyd on surface area, land-type-specific α values evaluated across seasons would be desirable to improve confidence in their use in global CTMs. As previously mentioned, Eugster and Hesterberg (1996) inferred a nocturnal non-stomatal median value for R_c(NO₂) of 700 s m⁻¹ (range 500–950 s m⁻¹) over managed grassland from soil-NO-corrected eddy covariance observations – similar to the median value of 750 s m⁻¹ simulated herein following Eq. (8) with α=1 (Fig. S8). Pilegaard et al. (1998) report nocturnal R_c(NO₂) of 771 ± 111 s m⁻¹ inferred from eddy covariance observations over a harvested wheat field (with re-growth) in southern Germany during mid-September. Given the high nocturnal NO₂ concentrations of 10–30 ppb at the location, soil NO most likely had a reduced effect on the resulting V_d(NO₂) compared to the large influence noted by Eugster and Hesterberg (1996), where nocturnal NO₂ concentrations were less than 10 ppb during periods of soil NO emission. The effect of soil NO on V_ex(NO₂) as a function of above-canopy NO₂ concentration can be seen in Fig. 2, where negligible influence is noted for concentrations above ∼ 5 ppb. Coe and Gallagher (1992) used eddy covariance to estimate a non-stomatal R_c(NO₂) of 548 s m⁻¹ over a heather moorland located in the southern Netherlands. Plake et al. (2015) find a maximum median nocturnal bulk R_c(NO₂) over a natural grassland site in Mainz, Germany, of 560 s m⁻¹ via the dynamic chamber approach. The nocturnal R_c(NO₂) values reported by Coe and Gallagher (1992) and Plake et al. (2015) are in between the r_hyd values computed here using α of 1 or 2. Assigning α=1 for low-roughness vegetative land types appears reasonable but may yield slight underestimates in nocturnal NO₂ uptake.

Reaction (R1) has been observed proceeding efficiently on ice surfaces, even at low temperatures (< 170 k) (Bang et al., 2015; Kim and Kang, 2010). Stocker et al. (1995) observed nocturnal deposition of NO₂ to a snow-covered grassland by eddy covariance in northern Colorado, reporting a median resistance to surface uptake of 740 s m⁻¹, from which we estimate an uptake coefficient to snow γ_snow of ∼ 1.6 × 10⁻⁵ (Tables 1 and S2) – similar to ${\mathit{\gamma }}_{\mathrm{g},{\mathrm{NO}}_{\mathrm{2}}}$ following Eq. (9) at 100 % RH. If NO₂ dry deposition persists into winter months at the levels observed for late fall in Fig. 3, this represents a significant depositional sink for wintertime NO₂ not currently represented in CTMs when both the lifetime and near-surface concentrations of NO_x are at a maximum.

3.3.4 Bottom-up estimates of nocturnal V_d(NO₂)

Simple estimates of bottom-up bulk-canopy V_d(NO₂) provide a useful sanity check on top-down eddy-covariance-inferred values and are a starting point for a mechanistic explanation of bulk-canopy deposition. Bottom-up estimates of nocturnal R_c(NO₂) for Harvard Forest were computed from parallel contributions of uptake to leaves, bark, and the forest floor:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{R}_{\mathrm{c}}\left({\mathrm{NO}}_{\mathrm{2}}\right)& =\left[\mathrm{1}/r_{\mathrm{c}}\mathrm{leaf}+\mathrm{1}/r_{\mathrm{c}}\mathrm{bark}\right.\\ & {\left.+\mathrm{1}/\left(r_{\mathrm{a}}+r_{\mathrm{c}}\mathrm{floor}\right)\right]}^{-\mathrm{1}},\end{array}\end{array}$$

where in-canopy aerodynamic resistance r_a was computed according to the Z03 scheme. Component canopy surface resistances in Eq. (10) were computed at an hourly resolution following Eq. (8) using the surface-specific NO₂ uptake coefficients from Table 1. Component surface area scale factors α are material dependent, varying according to the surface area used to normalize the deposition fluxes in corresponding measurement studies (Tables 1 and S2). Attention must be paid to the corresponding surface area over which a particular uptake coefficient is to be applied and to the understanding that application of an uptake coefficient for a complex surface where planar or simple estimates of geometric surface areas were used during measurement (i.e., forest floor, snow, or bark) assumes surface area equivalence in subsequent applications. As seen in Table 1, uptake of NO₂ to coniferous foliage under conditions of minimal stomatal aperture and RH < 90 % was found to be ∼ 2.6 times greater than to deciduous leaves on a projected area basis. In addition to the greater total-to-projected surface area of coniferous (∼ 2.7) compared to deciduous (∼ 2) leaves, we attribute the reduced non-stomatal uptake of NO₂ to deciduous leaves as a consequence of the hydrophobic adaxial (top) surface of deciduous leaves, which lack the elevated water vapor concentrations from stomata to support the thin water films for Reaction (R1) to proceed. We treat the leaves and bark of the forest canopy as wet following the dew flag from the Z03 scheme (Sect. 2.1.2). Under wet conditions, both top and bottom faces of deciduous leaves would be wetted, and we increase the α value used to scale uptake to deciduous leaves from LAI to 2LAI. The α value used to scale uptake to bark is computed as πSTAI, where STAI = 0.9 is the projected area index of tree branches (Horii et al., 2005). We compute forest floor surface resistances r_c floor as parallel contributions from uptake to snow and to the snow-free forest floor, weighted by the snow cover fraction following the Z03 scheme (Zhang et al., 2003a). As seen in Table 1, NO₂ uptake to wet bark is twice that to dry bark, and uptake to snow is approximately one-third that to the snow-free forest floor.

The top panel of Fig. 4 depicts monthly mean estimates of the nocturnal component canopy resistances used in Eq. (10) to compute bottom-up bulk-canopy R_c(NO₂). Also included for comparison are monthly mean estimates of top-down bulk-canopy r_hyd for α values of 1 and 2. The middle and bottom panels of Fig. 4 depict parameterized V_d(NO₂) following Eq. (2) using both bottom-up and top-down estimates of bulk-canopy surface resistance. Also included are the fractional contributions that leaf, bark, and forest floor surfaces make to total bottom-up NO₂ uptake. Due to the compensating seasonal contributions of r_c leaf and r_c floor to the total NO₂ deposition, bottom-up nocturnal V_d(NO₂) shows little seasonality, in accordance with eddy-covariance-inferred V_d(NO₂) from spring through fall (Fig. 3). As depicted in the middle panel of Fig. 4, bottom-up estimates of V_d(NO₂) computed using the uptake coefficients from Table 1 are greater than top-down optimized values across all seasons. Computing NO₂ uptake to dry and wet bark using the uptake coefficients from Table 1 of $\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{6}}$ and $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{5}}$, respectively, results in bark surfaces being the predominant nocturnal dry-depositional sink of NO₂ at Harvard Forest. Uptake coefficients to dry and wet bark in Table 1 are from Hanson et al. (1989), where chamber-measured NO₂ uptake was performed with Teflon end-capped cylindrical branch or trunk samples with diameters of ∼ 15 cm. It is expected that the calculated cylindrical geometric surface area used to normalize uptake to exposed bark samples in Hanson et al. (1989) is less than the surface area available for the reaction due to bark roughness and shape complexity for samples of this size. From Table 1, NO₂ uptake to smoother wood surfaces (i.e., wood board or plywood) at RH < 90 % is at least 3 times smaller than uptake to dry bark, which we speculate could result from the greater available surface area for the bark samples tested. Bark surface area for the trunk and branch samples examined in Hanson et al. (1989) may not be representative of the average bark surface area for the canopy at Harvard Forest. As depicted in the bottom panel of Fig. 4, reducing NO₂ uptake to bark by a factor of 2 (equivalent to a 2-fold increase in r_c bark) results in bottom-up V_d(NO₂) within the range of eddy-covariance-inferred monthly values corrected for soil NO (0.2–0.3 cm s⁻¹) and within 6 % of the top-down V_d(NO₂) for α=2 (parameterization P7) over a 12-month period. At this uptake level, bark is no longer the predominant sink for nocturnal dry deposition of NO₂, instead taking on a secondary role, where the predominant uptake alternates between canopy foliage during summer months and the forest floor during late fall, winter, and early spring. On an annual basis, nocturnal NO₂ uptake to forest surfaces is within 30 % for the reduced-bark-uptake case; specifically, the forest floor accounts for 37 % of uptake, bark 36 %, and foliage 27 %. Bottom-up modeling estimates of canopy-scale dry deposition of NO₂ would benefit from future chamber studies detailing NO₂ uptake to a variety of bark surfaces over a range of humidities and temperatures.

Figure 4 also depicts the nocturnal V_d(NO₂) computed using R_c(NO₂) parameterized following the Z03 scheme (Table S3). As with top-down and bottom-up estimates of V_d(NO₂), little seasonality is seen across monthly mean values. Over a 12-month period, V_d(NO₂) following Z03 is within 12 % and 17 % of top-down (α=2) and bottom-up (r_c bark×2) values, respectively. The Z03 scheme formulates non-stomatal uptake of NO₂ by scaling inverse surface resistances optimized for O₃ by 0.8 (Zhang et al., 2002a); herein, we explored the plausibility of non-stomatal NO₂ uptake being a result of Reaction (R1), with negligible contributions from other reactions. Although additional pathways for NO₂ uptake may exist in parallel to Reaction (R1) in some instances – particularly in leaf interiors (Sect. S6 in the Supplement) – non-stomatal NO₂ uptake to the natural land types examined and reviewed herein can plausibly be understood as occurring via Reaction (R1) when surface area effects are considered. Future work investigating the relative roles of Reaction (R1) and reduction reactions in NO₂ uptake to a variety of natural land type surfaces would aid mechanistic realism in future model developments.

Figure 4 (a) Computed estimates of monthly mean component canopy resistances for nocturnal NO₂ uptake to leaves, bark, and the forest floor at Harvard Forest and the resulting bottom-up bulk-canopy R_c following Eq. (10). Also included are top-down bulk-canopy resistances r_hyd following Eq. (8). (b) The resulting bottom-up and top-down canopy-scale nocturnal deposition velocities V_d(NO₂), including fractional contributions of leaf, bark, and forest floor surfaces to bottom-up V_d(NO₂). Also depicted is nocturnal V_d(NO₂) computed following the Z03 parameterization. (c) As in panel (b) but with deposition fractions and bottom-up V_d(NO₂) computed with uptake to bark reduced by a factor of 2. Monthly values of resistances, deposition velocities, and meteorological inputs are included in Table S3 in the Supplement. The surface-specific NO₂ uptake coefficients ${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$ used to calculate bottom-up R_c are included in Table 1.

3.4  Evaluation of parameterized V_d(NO_y) over Harvard Forest

Improved representation of V_d for NO_y component species HNO₃ and NO₂ via updates to molecular diffusivity (Sect. 3.2) and non-stomatal NO₂ uptake (Sect. 3.3.3), respectively, may be further evaluated at a broader scale through a full NO_y budget analysis. Simulated V_d(NO_y) from base and updated parameterizations is evaluated against observations from the HFEMS, considering NO_y component species NO₂, NO, HNO₃, and PAN. Of particular interest is the period from June to November 2000, when hourly observations of above-canopy HNO₃ concentration – a significant contributor to NO_y dry deposition at this location (Horii et al., 2005) – were added to the suit of measurements (Fig. S6). The top panel of Fig. 5 depicts the diel climatology of measured NO_y and component species over Harvard Forest alongside inferred NO_y computed as the sum of the component species. The middle panel of Fig. 5 depicts species-specific fractional contributions to measured NO_y. On average, NO_y inferred from the sum of measured component species is ∼ 76 % of measured total NO_y, with component contributions of 48 % (NO₂), 16 % (HNO₃), 8 % (PAN), and 4 % (NO).

To compare to V_d(NO_y) inferred from measured fluxes, we compute simulated deposition velocities V_d,sim(NO_y) from a linear combination of parameterized component deposition velocities V_d(x_i) weighted by species-specific concentration fractions (Michou et al., 2005; Wu et al., 2011):

$$\begin{array}{}\text{(11)}& V_{\mathrm{d},\mathrm{sim}}\left({\mathrm{NO}}_y\right)={\displaystyle \frac{\sum _i\left[x_i\right]V_{\mathrm{d}}\left(x_i\right)}{\sum _i\left[x_i\right]}}.\end{array}$$

Due to the large number of coincident hourly observations required for the comparison, only 19 coincident hourly values of V_d(NO_y) and V_d,sim(NO_y) exist across the entire dataset, which consists of over 2000 hourly measurement points. For this reason, a gap-filling method is used to estimate date- (d) and hour-specific (h) missing NO_y component concentrations ${\left[x_i\right]}_{d,h}^{\mathrm{infer}}$ for NO₂, NO, HNO₃, and PAN from measured NO_y:

$$\begin{array}{}\text{(12)}& {\left[x_i\right]}_{d,h}^{\mathrm{infer}}={\left({\displaystyle \frac{{\stackrel{\mathrm{‾}}{\left[x_i\right]}}_h^{\mathrm{meas}}}{{\stackrel{\mathrm{‾}}{\left[{\mathrm{NO}}_y\right]}}_h^{\mathrm{meas}}}}\right)}_{\mathrm{clim}}{\left[{\mathrm{NO}}_y\right]}_{d,h}^{\mathrm{meas}},\end{array}$$

where the diel climatologies of component fractions ${\stackrel{\mathrm{‾}}{\left[x_i\right]}}_h^{\mathrm{meas}}/{\stackrel{\mathrm{‾}}{\left[{\mathrm{NO}}_y\right]}}_h^{\mathrm{meas}}$ are depicted in the middle panel of Fig. 5. This method of gap filling was used by Wu et al. (2011) to evaluate simulated NO_y deposition velocities from the WRF-Chem and NOAH-GEM dry-deposition modules against this dataset. Here, we compute component fractions as the “ratio of smoothed means” rather than the “mean of ratios” to reduce the effect of outliers. The fraction of inferred species-specific hourly concentrations required for gap-filling over the study period is depicted in the bottom panel of Fig. 5 as the missing measurement fraction.

Diel profiles of the simulated V_d(x) used in Eq. (11) to compute simulated V_d(NO_y) are depicted in Fig. S4 and reviewed in Sect. S5 in the Supplement. Given the predominant contributions that NO₂ and HNO₃ make to the NO_y budget at Harvard Forest (Fig. 5), updates to the parameterization of V_d(NO₂) and V_d(HNO₃) discussed in previous sections will have a notable impact on simulated V_d(NO_y).

Figure 5 (a) Diel climatology of observed NO_y and component species NO₂, HNO₃, PAN, and NO measured above Harvard Forest (June–November 2000–2002). Solid lines depict a smoothing spline fit to hourly mean concentrations (dotted). Also shown is the sum of smoothed NO_y component species concentrations. (b) Fractional contributions of NO_y component species to observed NO_y. These values were computed as ratios of the smoothed diel mean concentrations in panel (a). Shown are individual NO_y species fractions (colored as above) and their sum (dashed black). (c) Fraction of hourly NO_y component species concentrations (colored as above) missing from the measured NO_y hourly time series spanning June–November 2000. As a gap-filling strategy for the calculated V_d(NO_y) hourly time series, missing values were inferred using Eq. (12).

3.4.1 Measurement–model comparison of V_d(NO_y)

Observations of hourly above-canopy eddy covariance fluxes of NO_y at the HFEMS are mostly downward (> 99 %; Fig. S9), regardless of adjustment for soil-emitted NO. Figure 6 depicts observation-inferred diel mean V_d(NO_y) alongside the simulated values from parameterizations P2, P4, and P8 – selected from Table 2 to highlight the dominant effects of diffusivity updates (P4) and surface resistance updates to NO₂ and PAN (P8) on simulated V_d(NO_y) referenced from the measurement height at Harvard Forest (P2). Inferred V_d(NO_y) was calculated from eddy-covariance-observed V_ex(NO_y) adjusted for estimates of soil-emitted NO, analogous to V_d(NO₂) in Eq. (6). As evident from Fig. 6, soil NO corrections to observed V_ex(NO_y) result in small increases to inferred V_d(NO_y) – 9 % for parameterizations P2 and P4, decreasing to 7 % for P8 due to an increase in the simulated canopy reduction factor (CRF) (Fig. S3) from updates to non-stomatal NO₂ uptake. By far the largest contributor to simulated V_d(NO_y) in Fig. 6 is HNO₃, contributing over 75 % to the dry deposition 24 h NO_y flux for parameterization P2 (Table 5). Despite nocturnal mean NO₂ concentrations on the order of 3 ppb (Fig. 5), NO₂ makes near-negligible contributions to simulated nocturnal V_d(NO_y) in parameterizations P2 and P4.

Figure 6 Simulated and observation-inferred diel mean NO_y dry-deposition velocities over Harvard Forest (June–November). In addition to observed V_ex(NO_y), inferred contributions to V_d(NO_y) from estimated soil NO emissions (Eq. 7) and maximum HONO emitted following Reaction (R1) are also depicted. Simulated V_d(NO_y) is depicted as the weighted sum (Eq. 11) of the contributing NO_y component species NO₂, HNO₃, and PAN for three simulation types: (P2) (base), (P4) (updated R_a and diffusivity D), and (P8) (updated R_c(NO₂) and R_c(PAN)).

Table 5 Inferred mean dry deposition fluxes of NO_y and component species over Harvard Forest.

^a See Table 2 for serial parameterization updates. Briefly, P2 is equivalent to the base P1 referenced from a Harvard Forest measurement height of 29 m. P3 has updates to aerodynamic resistance (R_a). P4 has updates to molecular diffusivity (D). P7 updates the resistance to surface uptake of NO₂ (R_c(NO₂)) to include Reaction (R1) using the surface area scale factor α=2 (Eq. 8). P8 updates the resistance to surface uptake of PAN (R_c(PAN)). ^b 24 h mean (±1σ) dry deposition fluxes inferred from the product of simulated V_d(x) and gap-filled measured concentrations over June–November 2000. ^c Inferred NO_y flux from the sum of inferred component (HNO₃, NO₂, and PAN) fluxes.

Updates to the parameterization of molecular diffusivity (P4) results in large reductions in simulated V_d(NO_y) (−30 % during the day, −28 % at night) and consequently in a 24 h depositional flux (26 % reduction, Table 5) due to large reductions in simulated V_d(HNO₃) (Fig. S4). This update exposes a morning peak in observation-inferred V_d(NO_y) in Fig. 6, which simulated values fail to capture. Given the heavy reliance on the climatological diel profile of HNO₃ (Fig. 5) due to the paucity of hourly HNO₃ observations (Figs. 5 and S6) and the potential influence of unmeasured rapidly depositing NO_y species (Horii et al., 2005), further discussion of this model–measurement discrepancy is beyond the scope of this work; however, we note that it overlaps with the period of mixed layer growth.

As seen for parameterization P8 in Fig. 6, updates to non-stomatal deposition of NO₂ and PAN result in noticeable increases in simulated V_d(NO_y) (8 % during the day, 30 % at night), with updates to NO₂ deposition having a greater relative impact at night (25 % increase in V_d(NO_y)) than during the day (5 % increase in V_d(NO_y)). Parameterization P8 updates are associated with an 18 % increase in the inferred 24 h NO_y dry deposition flux due to large increases in inferred dry deposition of NO₂ (56 % increase) and PAN (85 % increase) (Table 5).

Also depicted in Fig. 6 is the effect on observation-inferred V_d(NO_y) from a maximum estimate of HONO emitted from the heterogeneous hydrolysis of NO₂ on deposition surfaces following Reaction (R1). Assuming unity for HONO surface emission to the gas phase and subsequent ventilation from the canopy produces a 4 % increase in inferred V_d(NO_y). Uncertainties exist regarding the nature of the dynamic equilibrium that establishes between evolved and adsorbed/dissolved HONO (Collins et al., 2018; Harrison et al., 1996; Lee, 2012; Spicer et al., 1993; Wojtal et al., 2011) and subsequent implications of a nocturnal reservoir of deposited HONO as a daytime source of HONO to the atmospheric surface layer (He et al., 2006; Ren et al., 2020; VandenBoer et al., 2014, 2015; Wentworth et al., 2016). Lee (2012) monitored near-continuous above-canopy NO₂ and HONO concentrations and eddy covariance fluxes at the HFEMS during 2011, finding nocturnal enhancements in HONO concomitant with NO₂. However, neither upward nor downward fluxes of HONO were observed, suggesting the establishment of dynamic equilibrium between HONO emission and deposition at Harvard Forest, where perturbation fluxes were below detection limits. Measurements of HONO and NO₂ fluxes over grassland and sugar beet surfaces have highlighted the bidirectional nature of HONO exchange; HONO emission was found to dominate the bidirectional flux under elevated NO₂ concentrations (> 10 ppb for the land types studied), while deposition was noted to be dominant at lower ambient NO₂ concentrations (Harrison et al., 1996; Harrison and Kitto, 1994). The diel flux behavior of HONO is likely multifactorial, depending on land type, meteorology, and trace gas and particulate concentrations (Pusede et al., 2015; VandenBoer et al., 2015). The absence of significant fluxes of HONO over Harvard Forest, as noted by Lee (2012) despite observed downward nocturnal fluxes of NO₂ and nocturnal enhancement of HONO, does not exclude Reaction (R1) as a source of HONO. Rather, this indicates the importance of deposition and re-emission processes from canopy surfaces, which may dominate HONO behavior at rural forest sites (Ren et al., 2011, 2020; Sörgel et al., 2011; Zhou et al., 2002) in the absence of strong pulses of ambient NO₂ perturbing the dynamic equilibrium between adsorbed/dissolved and gas-phase HONO – as is routinely observed in laboratory studies (Finlayson-Pitts et al., 2003; Spicer et al., 1993).

Neglecting the possible effect of emitted HONO on V_d(NO_y), given the findings of Lee (2012), we find a moderate underestimate of 20 % in simulated V_d(NO_y) (24 h) from parameterization P8, with similar daytime and nighttime biases of −20 % and −19 %, respectively. As previously discussed in Sect. 3.4 and depicted in Fig. 5b, the NO_y concentration budget at the HFEMS is closed to 76 % on average from observations of NO_x, HNO₃, and PAN. Horii et al. (2005) provide evidence of a rapidly depositing unidentified NO_y species at this site, especially under southwesterly flow, and speculate that the unidentified NO_y species are organic nitrates.

4 Assumptions of nocturnal stomatal closure and non-stomatal NO₂ uptake via Reaction (R1)

Despite long-standing uncertainty regarding nocturnal stomatal behavior (Caird et al., 2007; Costa et al., 2015; Dawson et al., 2007), it is generally assumed that at night in response to elevated guard cell CO₂ concentration, the stomata of C₃ and C₄ plants are nearly closed (Nobel, 2009), thereby restricting trace gas exchange. In both the widely used W89 and Z03 dry-deposition schemes, stomata are assumed to be closed at night, and therefore all nocturnal deposition is parameterized through non-stomatal pathways. Although chamber studies consistently find that stomatal control dominates daytime foliar uptake of NO₂, uncertainty remains regarding the importance of non-stomatal pathways at night. Some studies note negligible non-stomatal contributions to NO₂ deposition, attributing nocturnal uptake to partially open stomata (Chaparro-Suarez et al., 2011; Delaria et al., 2020; Geßler et al., 2000; Rondón et al., 1993; Sparks et al., 2001), while others find non-stomatal contributions to be non-negligible (Geßler et al., 2002; Hanson et al., 1989; Thoene et al., 1996; Wang et al., 2020; Weber and Renenberg, 1996). This amalgam of contradictory results warrants further consideration and is discussed in the context of a literature review in Sect. S6 in the Supplement.

Notwithstanding nocturnal stomatal behavior, nocturnal ambient humidity is often elevated at locations with lush vegetation, as it is during summer months at Harvard Forest (monthly median RH 92 %–95 %; Table S3), whereby growth of aqueous films has been observed to occlude stomatal pores (Grantz et al., 2018); heterogeneous hydrolysis of NO₂ would be expected to proceed on foliar as well as other surfaces under these conditions. Furthermore, nocturnal uptake of NO₂ through a deposition process involving HONO production has been implicated in the field on several occasions (Harrison and Kitto, 1994; Ren et al., 2020; Stutz et al., 2002), including with RH dependence (Stutz et al., 2004; VandenBoer et al., 2013).

Inferred values of ${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$ in Tables 1 and S2 fall within the range expected for uptake due to NO₂ heterogeneous hydrolysis, i.e., 10⁻⁶ to 10⁻⁵ (Sect. S3). For foliar surfaces under dark conditions, “total leaf area” normalized uptake was not observed to exceed that of the planar surface of distilled water, supporting the possibility that Reaction (R1) is the predominant mechanism driving uptake to both the non-foliar and foliar surfaces presented. Three features stand out from the tabulated values of ${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$ in Tables 1 and S2. First, there is a dependence on surface moisture – dependence on RH or surface wetness is seen for uptake to wood board, concrete, and tree bark. The RH or wetness dependence for NO₂ uptake to wood board and bark is similar to the factor of 2 increase in ${\mathit{\gamma }}_{\mathrm{g},{\mathrm{NO}}_{\mathrm{2}}}$ from Eq. (9) between RH 50 % and 100 %. Second, the surface area available for heterogeneous reaction has a direct influence on the resulting surface-specific uptake, as expected for a collision-limited heterogeneous process such as Reaction (R1). Surfaces with complex and undetermined microscopic surface areas (i.e., bark, coarse concrete, forest floors, and snow) exhibit higher $V_{\mathrm{d}}^{\mathrm{surf}}$ and resulting ${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$ – a factor of 2 to 30 greater than surfaces normalized by accurate predictions of available surface area (i.e., bulk water and foliar). This increased uptake to convoluted surfaces could be an indirect measure of the total microscopic surface area available for reaction, thus highlighting the utility of using field-derived uptake coefficients to parameterize surface uptake in dry-deposition models. Third, there is a distinctive feature between NO₂ uptake to coniferous and deciduous leaves when stomatal aperture is at a minimum under dark conditions. When uptake is normalized by projected leaf area, inferred ${\mathit{\gamma }}_{{\mathrm{NO}}_{\mathrm{2}}}$ to coniferous species is on average 2.6 times greater than to deciduous species – a factor similar to the ratio of total-to-projected needle leaf area (Riederer et al., 1988). Since coniferous needles have stomata distributed across the entire leaf surface, whereas most deciduous leaves have stomata located on the lower (abaxial) leaf surface and a thicker hydrophobic wax cuticle on the upper (adaxial) leaf surface, the ∼ 2.6-fold greater uptake to coniferous species (when normalized by projected leaf area) may reflect the absence of thin water films on the adaxial surface of deciduous leaves in chamber studies from Table 1. These inferences are consistent with the work of Sumner et al. (2004), where similar rates of NO₂ heterogeneous hydrolysis across a variety of hydrophilic and hydrophobic surfaces were understood in the context of the available surface areas supporting thin water films. However, given the intra- and interspecies variability in nocturnal $V_{\mathrm{d}}^{\mathrm{surf}}$ from leaf-level studies across a range of environmental conditions (i.e, Breuninger et al., 2013; Delaria et al., 2018, 2020; Rondón et al., 1993; Wang et al., 2020), further investigation into this simple generalization is warranted.

5 Conclusions

A stand-alone version of the gaseous dry-deposition algorithm from GEOS-Chem, implemented to run in single-point mode, enabled a detailed and more direct evaluation of various branches of the algorithm against eddy-covariance-inferred deposition velocities over two North American temperate forest ecosystems. Observations of deposition velocities for species that deposit under dynamical control (i.e., nominally small surface uptake resistance R_c) facilitated the identification of a large high bias in computed molecular diffusivities. Correction of this bias using Fuller's method to calculate diffusivities in the absence of measured values resulted in improved simulation of V_d(HNO₃). Consequently, this exposed a morning peak in observation-inferred V_d(NO_y) that simulated values failed to capture. Site-specific roughness length and reference height were found to be important constraints on the calculation of aerodynamic resistance for rapidly depositing species, whereas correction for the influence of the roughness sublayer was found to be of minor importance at the measurement heights involved and to have a negligible effect at the dry-deposition reference height used in GEOS-Chem, in agreement with previous work (Simpson et al., 1998).

A large low bias (−80 %) in simulated nocturnal V_d(NO₂) against eddy-covariance-inferred values spanning many months over Harvard Forest was identified in accordance with previous work (Horii, 2002; Horii et al., 2004). We addressed this low bias by representing NO₂ heterogeneous hydrolysis – a well-known surface reaction (Finlayson-Pitts et al., 2003) shown to be of importance in the field (VandenBoer et al., 2013) but to our knowledge still unrepresented in dry-deposition parameterizations – in the calculation of non-stomatal surface resistance. A literature review of surface-specific NO₂ deposition velocities to both non-foliar and nocturnal foliar surfaces highlights the importance of considering microscopic surface area for heterogeneous reaction and enabled estimates of bottom-up V_d(NO₂) for Harvard Forest, which agree well with top-down estimates from eddy-covariance-inferred values when uptake to bark was reduced by a factor of 2 – indicating a need for further laboratory study of uptake to bark for representative samples. Consideration of the effect of soil NO emission on eddy-covariance-inferred V_d(NO₂) at Harvard Forest was found to be important, as was representative canopy surface area when applying NO₂ uptake coefficients from field and laboratory observations.

We persist with the assumption that nocturnal uptake of NO₂ follows non-stomatal pathways, as is currently the case in the dry-deposition schemes widely used in atmospheric CTMs. Meanwhile, the nocturnal behavior of stomata remains an active area of research. Confounding processes such as the hydraulic activation of stomata (HAS) and condensation at elevated RH complicate the inference of stomatal conductance to trace gases from observations of water vapor flux, especially in well-mixed chambers under dark conditions. It would be helpful for future enclosure studies of NO₂ uptake to consider the effects of heterogeneous hydrolysis of NO₂ on foliar and non-foliar surfaces, as well as potential biases in estimates of stomatal conductance resulting from possible HAS.

Although we find that the Z03 dry-deposition scheme adequately captures the magnitude of nocturnal V_d(NO₂) over Harvard Forest, formulating uptake relative to V_d(O₃) neglects contributions from Reaction (R1) and therefore may be a misrepresentation for NO₂, with implications for surface HONO production in models. We recommend consideration of a non-stomatal dry-deposition scheme for NO₂ in atmospheric models that considers the effect of NO₂ hydrolysis on surface resistances. This represents a significant depositional sink for NO₂ under conditions when both the lifetime and the near-surface concentration of NO₂ may be elevated, i.e., nocturnal and urban wintertime conditions. We present two approaches that result in general agreement for a mature temperate forest ecosystem. The simplest approach is to represent non-stomatal resistance to NO₂ uptake as r_hyd following Eq. (8), with a surface area scale factor α=2 for high-surface-area land types such as urban and forest, and α=1 for remaining land types. We caution that the value α=2 is not a fit to observations at Harvard Forest; instead, it is likely an upper estimate from a sensitivity study. Long-term field studies quantifying atmosphere–surface exchange across a variety of land types and seasons would facilitate further development of this unique dry-deposition pathway.

Although a main focus of this work was on the effect that NO₂ hydrolysis has on V_d(NO₂), there is much interest in accurate HONO simulation, given that the near-surface nocturnal buildup of HONO results in an early-morning burst of OH and NO radicals as HONO photolyzes (Finlayson-Pitts, 2009; Ren et al., 2020), initiating photochemistry prior to other HO_x precursors (Platt et al., 1980). Future work will present the implementation of updates to V_d(NO₂) developed in this paper within the GEOS-Chem CTM, including analysis of impacts on simulated concentration fields (i.e., NO₂, HONO, OH, and O₃) and deposition budgets.