Sensitivity of simulated ammonia fluxes in Rocky Mountain National Park to measurement time resolution and meteorological inputs

Gaseous ammonia (NH₃) is an important precursor for secondary aerosol formation and contributes to reactive nitrogen deposition. NH₃ dry deposition is poorly quantified due to the complex bidirectional nature of NH₃ atmosphere-surface exchange and lack of high time-resolution in situ NH₃ concentration and meteorological measurements. To better quantify NH₃ dry deposition, measurements of NH₃ were made above a subalpine forest canopy in Rocky Mountain National Park (RMNP) and used with in situ micrometeorology to simulate bidirectional fluxes. NH₃ dry deposition was largest during the summer, with 47 % of annual net NH₃ dry deposition occurring in June, July, and August. Because in situ, high time-resolution concentration and meteorological data are often unavailable, the impacts on estimated deposition from utilizing more commonly available biweekly NH₃ measurements and ERA5 meteorology were evaluated. Fluxes simulated with biweekly NH₃ concentrations, commonly available from NH₃ monitoring networks, underestimated NH₃ dry deposition by 45 %. These fluxes were strongly correlated with 30 min fluxes integrated to a biweekly basis (R² = 0.88), indicating that a correction factor could be applied to mitigate the observed bias. Application of an average NH₃ diel concentration pattern to the biweekly NH₃ concentration data removed the observed low bias. Annual NH₃ dry deposition from fluxes simulated with reanalysis meteorological inputs exceeded simulations using in situ meteorology measurements by a factor of 2.

1 Introduction

Gaseous ammonia (NH₃) is an important atmospheric constituent, with effects on atmospheric chemistry and the nitrogen cycle. Atmospheric deposition of reactive nitrogen (N_r) is linked to nitrogen oxides (NO_x) and NH₃ emissions. Emissions of NO_x and NH₃ have many potential fates, including chemical transformation, dry deposition, particle formation, and wet deposition. Anthropogenic emissions of NO_x and NH₃ are produced predominantly by combustion and agriculture, respectively (Walker et al., 2019a), although there are also NH₃ emissions from traffic, wastewater treatment, and wildfires (Walker et al., 2019b; Tomsche et al., 2023). Due to increased food demand and industrialization, anthropogenic N_r has been increasing annually (Galloway et al., 2008; Kanakidou et al., 2016). Excess reactive nitrogen deposition has well-documented adverse effects on ecosystem health, including eutrophication, soil acidification, decreased biodiversity, and increased N in freshwater bodies (Bobbink, 1991; Baron, 2006; Holtgrieve et al., 2011; Boot et al., 2016; Zhan et al., 2017; Pan et al., 2021).

As a result of effective NO_x emission controls, the balance of N_r wet deposition across the US has shifted from oxidized N-dominated to reduced N-dominated, and dry deposition of NH₃ at times dominates total N_r deposition (Li et al., 2016; Walker et al., 2019a; Driscoll et al., 2024). The National Emission Inventory (NEI) indicates that US NO_x emissions were reduced by 46 % between 2013 and 2023, while NH₃ emissions increased by 13 % (U.S. EPA, 2023).

Critical loads, deposition levels below which harmful effects are not expected to occur, have been estimated for many ecosystems (e.g., Bowman et al., 2012; Schwede and Lear, 2014). In Rocky Mountain National Park (RMNP), a critical load of 1.5 kg N ha⁻¹ yr⁻¹, based on wet deposition of NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$, has been established to avoid adverse effects on the ecosystem (Baron, 2006). The pre-industrial nitrogen load has been estimated at 0.2 kg N ha⁻¹ yr⁻¹, while the current wet deposition rate is as high as 3.65 kg N ha⁻¹ yr⁻¹, approximately 15× the natural background and significantly higher than the critical load (Burns, 2003; CDPHE, 2007; Benedict et al., 2013a). Although the RMNP N_r critical load only considers wet deposition of NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$, dry deposition can also contribute significantly to total N_r deposition. NH₃ dry deposition in RMNP was estimated to be the third largest contributor to total N_r deposition, accounting for 18 % of N_r deposition from November 2008 to November 2009 (Benedict et al., 2013a).

NH₃ dry deposition, however, remains a highly uncertain component of N_r deposition, and fluxes are rarely measured (Walker et al., 2019b). Previous studies in RMNP have estimated NH₃ dry deposition using unidirectional inferential models, where the NH₃ deposition velocity (V_d) was approximated as 70 % of the HNO₃ deposition velocity (Beem et al., 2010; Benedict et al., 2013a, b) and NH₃ emission from the surface was ignored. In reality, NH₃ exchange between the atmosphere and surface is bidirectional, including deposition to and emission from the surface (Sutton et al., 1995). Several models have been developed to simulate the bidirectional exchange of NH₃ with the surface (Massad et al., 2010; Zhang et al., 2010; Pleim et al., 2013). Key model inputs include micrometeorology, soil and vegetation parameters, and atmospheric concentrations. In practice, fluxes can change quickly and even reverse direction with changing environmental conditions. Gaseous NH₃ is challenging and expensive to measure at high time resolution; lower-cost weekly or biweekly passive diffusion-based sampler measurements are more commonly utilized for long-term monitoring (Li et al., 2016; Schiferl et al., 2016; Butler et al., 2016; Hu et al., 2021). Previous efforts have used these low-cost measurements to quantify NH₃ dry deposition (Walker et al., 2008; Shen et al., 2016; Tanner et al., 2022). Detailed, high time-resolution meteorological observations at the location of interest are also desired when estimating dry deposition. Due to the frequent unavailability of such data, reanalysis meteorological data is often used as a substitute (Wichink Kruit et al., 2012; Schrader et al., 2018).

Schrader et al. (2018) investigated the impact of low time-resolution NH₃ concentrations on modeled fluxes. They found that using monthly NH₃ concentrations underestimates total NH₃ dry deposition. However, due to a linear relationship between simulations using monthly NH₃ concentrations and those using hourly NH₃ concentrations, they were able to generate a site-specific correction to compensate for the use of low time-resolution concentration data. Simulations were done using a simplified parameterization of the bidirectional exchange model described in Massad et al. (2010), and the NH₃ concentrations were simulated using the LOTOS-EUROS model (Hendriks et al., 2016).

Understanding and managing these biases could unveil opportunities to estimate NH₃ deposition when high time-resolution, in situ concentration, and meteorological observations are unavailable. Using high time-resolution NH₃ concentration measurements, we provide the first estimate of NH₃ annual dry deposition to an RMNP forest canopy using a bidirectional exchange model driven by high time-resolution NH₃ concentration data and in situ micro-meteorological measurements. We use in situ data collected in RMNP to determine if site-specific correction factors suggested by Schrader et al. (2018) apply to real-world observations and whether correction factors can be employed to reduce biases associated with NH₃ simulations using lower-cost, low time-resolution NH₃ measurements such as those available from the U.S. Ammonia Monitoring network (AMoN) (Puchalski et al., 2011). We also tested if an average NH₃ diel pattern could be applied to reduce these biases and, if so, the length of measurements necessary to adequately describe the diel pattern. Finally, we examine biases introduced by substituting reanalysis meteorological data for high time-resolution in situ measurements.

2 Data and methods

2.1 Site location

Study observations were collected in RMNP in northern Colorado. The park, established to preserve the natural landscape, including montane, subalpine, and alpine ecosystems, is predominantly above 3000 m, where ecosystems developed under nutrient-deprived conditions and are therefore sensitive to excess inputs of nitrogen. Nitrogen deposition has been a historical problem in RMNP, with diatom changes documented starting in the 1950s and more recent effects, including eutrophication and changes to plant species (Baron et al., 2000; Korb and Ranker, 2001; Baron, 2006).

The area east of RMNP (Fig. 1) includes a large urban corridor and extensive agricultural activity in the plains. The Front Range urban corridor, spanning from Denver to Fort Collins, is a major source of nitrogen oxide emissions (Benedict et al., 2013b). The northeast plains of Colorado are predominantly agricultural and include major sources of NH₃ emissions from both animal feeding operations and crop production. The spatial pattern seen for feedlots is broadly consistent with the spatial distribution of other agricultural activities. Pan et al. (2021) found that 40 % of summertime dry deposition of NH₃ in RMNP was associated with transport from agricultural regions to the east.

Figure 1 A map of the study region. Animal units are shown as the number of permitted animals as of 2017, scaled by an animal unit factor relative to the species. Elevation data is from the US Geological Survey Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010) at 7.5 arcsec spatial resolution, or 225 m (available at: https://earthexplorer.usgs.gov/, last access: 27 February 2024).

Data was collected at two adjacent locations for this study, both near the base of Longs Peak in Rocky Mountain National Park: a National Ecological Observatory Network (NEON) tower site (40.275903, −105.54596) and a nearby National Park Service shelter (∼ 500 m north of the NEON tower), from September 2021 through August 2022. The study location, denoted with a star in Fig. 1, is 2750 m above sea level. The tower is surrounded by lower montane forest, comprised of predominantly evergreen needleleaf species, including ponderosa pine, juniper, and Douglas fir. There are also groves of quaking aspen located in the region. Meteorological transport to the site is generally bimodal. Prevailing downslope transport from the northwest occurs generally overnight and during the cooler months, when ammonia concentrations are typically low. The mountain-plains circulation generates daytime upslope transport, bringing air masses from the plains east of the park up into RMNP. This pattern strengthens during warmer seasons. Periods of synoptically forced sustained upslope transport are also common, especially during spring and autumn (Gebhart et al., 2011). Downslope and upslope transport patterns are not due west and east at the study site because of channelling by local topography.

At RMNP, a diel pattern in ambient NH₃ concentrations has commonly been observed in past measurements. This pattern is primarily driven by changes in transport patterns that carry NH₃ emissions to the park (Benedict et al., 2013b; Juncosa Calahorrano et al., 2024) and, sometimes, modified by changes in the atmosphere-surface exchange of NH₃, especially during NH₃ uptake and emission from dew formation and evaporation (Wentworth et al., 2016).

2.2 Micrometeorological measurements

2.2.1 In situ micrometeorology

Meteorological and soil data were accessed from the RMNP NEON flux tower. The mean canopy height in the area surrounding the tower is 19 m. Temperature (mean = 6 °C), relative humidity (mean = 40 %), and annual days of precipitation are highly variable at the site due to its high elevation. Mean values were calculated from September 2021 to September 2022. Snowfall typically occurs between October and May. The seasonal mean temperatures (relative humidities) are as follows: winter (December, January, and February) mean is −3 °C (30 %), spring (March, April, and May) mean is 2 °C (44 %), the summer (June, July, and August) mean is 15 °C (49 %), and the fall (September, October, and November) mean is 8 °C (37 %). Precipitation is measured at 1 min resolution by a Belfort AEPG II 600M weighing gauge. Precipitation events were defined as periods of rainfall separated by at least one hour without precipitation. During our study period, there were 27 precipitation events in the winter, 62 in the spring, 63 in the summer, and 26 in the fall.

Meteorological data accessed from the NEON site includes wind vectors, friction velocity, Obukhov length, soil temperature, short wave radiation, relative humidity, air density, air pressure, and air temperature above the tree canopy. The meteorological observations used from the NEON tower are 30 min mean values. Direct measurements of wind vectors, air temperature, short wave radiation, relative humidity, air density, and air pressure were used from the tower-top measurements (25 m a.g.l., meter above ground level). 3D wind vectors were measured at 20 Hz using the CSAT-3 sonic anemometer (Campbell Scientific Inc., Logan, Utah, USA). Soil temperature was taken as the average across 5 collection sites within 200 m of the flux tower. Leaf area index (LAI) is estimated at the site using remotely sensed data. The square kilometer of leaf area index values surrounding the tower site is shown in Fig. S5 in the Supplement. A mean value of 0.8 was estimated using the landscape surrounding the site. The sensitivity of simulated NH₃ fluxes to LAI can also be found in Sect. S5 of the Supplement. Additional information about each of the reported NEON datasets can be found in the Site Management and Event Reporting documentation (available at: https://doi.org/10.48443/9p2t-hj77, NEON, 2023).

NEON meteorological data contained gaps due to power outages and scheduled instrument maintenance. Across the year of data, the gaps comprised 5.8 % of the data (1021 data points). To quantify the annual deposition of NH₃ in RMNP, these gaps were filled using the average diel pattern of fluxes during the current biweekly NH₃ sampling period.

2.2.2 Reanalysis meteorology data

Detailed meteorological and soil data are not available at many locations where NH₃ dry deposition is of interest. Reanalysis data, which combine short-range weather forecasts with assimilated observations, are a common source of meteorological data that can be used in the absence of local observations. To probe the impact of using reanalysis data in place of in situ observations, a set of bidirectional flux simulations was conducted using ERA5 hourly reanalysis data (Hersbach et al., 2020). ERA5 hourly reanalysis data has a spatial resolution of 0.25°, or approximately 31 km. The parameters used from the ERA5 data are as follows: air temperature, air pressure, dewpoint temperature, turbulent surface stress, moisture flux, sensible heat flux, friction velocity, standard deviation of filtered subgrid orography, solar radiation, and soil temperature. Obukhov length (L) is not given in the ERA5 dataset and was calculated using Eq. (1) following Stull (1988), shown below. Obukhov Length is the characteristic length scale of turbulence in the atmospheric boundary layer and is calculated from ERA5 data using instantaneous surface sensible heat and moisture fluxes based on the suggested calculation from the European Centre for Medium-Range Weather Forecasts (Guisti, 2024).

$$\begin{array}{}\text{(1)}& L={\displaystyle \frac{-\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{\mathrm{v}}^{\prime }}{u}_{\ast }^{\mathrm{3}}}{kg\stackrel{\mathrm{‾}}{{\left(w^{\prime }{\mathit{\theta }}_{\mathrm{v}}^{\prime }\right)}_s}}},\end{array}$$

where k is the von Karman constant, g is gravitational acceleration, $\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{\mathrm{v}}^{\prime }}$ is the mean virtual temperature near the surface, $\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}_{\mathrm{v}}^{\prime }}$ is the surface flux of virtual potential temperature, and u_∗ is the friction velocity.

2.3 NH₃ data

2.3.1 Biweekly NH₃ measurements

Biweekly NH₃ ambient air concentration was measured using Radiello passive diffusion samplers purchased from Sigma Aldrich. The Radiello sampling system includes a diffusive body (part number: RAD1201) and an adsorbing cartridge (part number: RAD168) coated with phosphoric acid. NH₃(g) diffuses across the exterior diffusive body and is collected on the adsorbing cartridge as ammonium (NH${}_{\mathrm{4}}^{+}$) over two weeks. Collected ammonia (as NH${}_{\mathrm{4}}^{+}$) is extracted from the cartridge into deionized water and analyzed on a cation IC using a 20 mM methanesulfonic acid eluent (0.5 mL min⁻¹) on a Dionex CS12A ion exchange column with a CSRS ULTRA II suppressor and Dionex conductivity detector (Li et al., 2016). NH₃ passive samples were collected in duplicate ($\mathit{\sigma }=\pm \mathrm{0.25}$ µg m⁻³) on top of the NEON tower (25.35 m a.g.l.). Across the study period, there were 27 sampling periods. Due to site access issues, some samples had durations longer than 2 weeks. To create a consistent dataset, all data were aggregated to a 2-week average. In the case where two samples overlapped during a 2-week period, they were combined using a weighted average. One sample was below the detection limit and was removed from this analysis. Passive NH₃ sampling methods have been shown to have a low bias when compared with other sampling methods, including annular denuders and Picarro Cavity Ringdown spectroscopy methods (Puchalski et al., 2011; Pan et al., 2020).

2.3.2 High temporal resolution NH₃ measurements

NH₃(g) air concentration was also measured using an ion mobility spectrometer (IMS). Ion mobility spectroscopy separates ionized molecules based on their mobility through a carrier gas, under the influence of an electric field. The instrument used was the AirSentry II Point-of-Use IMS (Particle Measuring Systems, Niwot, CO). The instrument was in the National Park Service (NPS) shelter (located at 40.278129, −105.545635), 500 m north of the NEON site, with an inlet located approximately 2 m above natural grassland. The sampling inlet was 0.635 cm Teflon tubing, heated to 40 °C to reduce NH₃ loss to the sampling tube. Inlet length was kept as short as possible to further prevent NH₃ loss. Particles were removed by a fiber filter at the tip of the inlet. Due to the high altitude of the site location, the instrument was zeroed to account for pressure differences upon installation. Multi-point calibrations were conducted at the beginning and end of sampling. Calibration was confirmed using a known concentration ammonia gas sample split between the instrument and a phosphoric acid-coated annular denuder, where the NH₃ collected by the denuder is extracted into deionized water and analyzed using ion chromatography. Zero measurements were made periodically by overflowing the inlet with ultra-high purity clean air. The AirSentry samples at a 30 s frequency. During the study, the AirSentry collected 919 000 data points. The limit of detection for 30 s measurements 70 pptv. For this data analysis, NH₃ concentration data was averaged to 30 min mean values. Averaging data points increases the signal-to-noise ratio. We approximate that the signal-to-noise ratio increases proportionally to the square root of the number of samples (n=60) (Dempster, 2001). In this case, the signal-to-noise ratio increases by a factor of 7.7, reducing the limit of detection to 9 pptv for 30 min mean NH₃ concentrations. Across the year of data collection, 101 points fell below the detection limit.

2.3.3 NH₃ data preparation

To investigate the effect of NH₃(g) sampling time-resolution on simulated fluxes, bidirectional fluxes were simulated with concentration data at: (i) 30 min resolution (30 min NH₃), (ii) with the 2-week integrated passive NH₃ (Biweekly Passive NH₃), and (iii) lastly with an average diel profile applied to each day within the 2-week passive period (Average Diel Pattern NH₃). The three NH₃ data products are shown in Fig. 2.

Figure 2 Three NH₃ concentration data sets are shown for the entire study period. The two-week average across each concentration data product is the same.

The 30 min NH₃ concentration data is generated using a combination of data from the AirSentry NH₃ located at the NPS shelter and passive NH₃ samples collected on the NEON tower. Data gaps, due to power outages and regular maintenance, were filled using the average diel pattern across the year of data collection. Data gaps accounted for about 3000 out of more than 900 000 points across the study period. To generate a 30 min NH₃ data set above the tree canopy, the data was divided into biweekly periods that match the passive NH₃ collection periods. The average concentration from the AirSentry across each period was then scaled to match the biweekly passive NH₃ concentration. The 101 30 min average NH₃ concentration values below the AirSentry detection limit, representing 0.5 % of the total measurement period, were assumed to represent a random distribution below the detection limit and retained for post-process scaling from the passive observations. This preserves the temporal variability of NH₃ concentrations while ensuring that the average air concentration across the sampling period is consistent with the passive NH₃ measurements atop the NEON tower, which can differ from those above the adjacent grassland where the Air Sentry measurements are made.

The biweekly passive NH₃ with diel profile applied is generated using the annual average diel pattern of NH₃ from the AirSentry data. To determine if there are systematic differences between the NH₃ diel pattern at the two sites, raw and scaled AirSentry concentrations were compared to 4 and 6 h University Research Glassware denuder measurements taken on the NEON tower. The NH₃ concentrations were well correlated between sites. This comparison is shown in Fig. S1. Each day of the biweekly passive period is assigned the average diel pattern, then the biweekly mean is scaled to match the biweekly passive concentration. This dataset was generated to investigate if the inclusion of a simple diel profile was sufficient to correct for the bias in bidirectional fluxes created by using low time-resolution NH₃ concentrations, as shown by Schrader et al. (2018).

These three concentration data sets will be used for bidirectional flux simulations of NH₃. For the rest of this work, the three NH₃ data sets will be referred to using the following nomenclature.

•  

  30 min NH₃: NH₃ concentration data at 30 min resolution

•  

  Biweekly NH₃: Biweekly Passive NH₃ concentration data

•  

  Average Diel Pattern NH₃: Passive NH₃ concentration scaled using an average diel profile from the 30 min NH₃ dataset

2.4 Additional measurements

2.4.1 Wet deposition data

Weekly precipitation wet deposition data was obtained from the National Trends Network (NTN) (National Atmospheric Deposition Program, 2022) site at Beaver Meadows in RMNP (“CO19”: located at 40.3639, −105.5810). The Beaver Meadows site location, at 2477 m elevation and located approximately 10 km north of the CASTNET site, is shown in Fig. 1.

2.4.2 Additional gas and particle measurements

Additional air concentration data was obtained from the U.S. EPA Clean Air Status and Trends Network (CASTNET) site at the NPS shelter (“ROM206”: located at 40.278129, −105.545635). Weekly filter pack concentrations of nitric acid (HNO₃) and sulfur dioxide (SO₂) were used to calculate the acid ratio (Eq. 10) in the bidirectional exchange simulations of NH₃ (U.S. EPA, 2024a).

Weekly dry deposition of HNO₃, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ was estimated by CASTNET (U.S. EPA, 2024b) using the weekly filter pack concentrations and historical values of deposition velocity (V_d) from the U.S. EPA Multi-Layer Model (MLM) (Meyers et al., 1998). The generation of deposition velocities was discontinued in 2019. Bowker et al. (2011) found that using historical values of V_d from the U.S. EPA Multi-Layer Model did not significantly bias the annual mean of deposition.

One approach to estimating NH₃ deposition is to estimate the V_d as a fixed fraction (70 %) of the V_d of HNO₃. This approach has been historically used to estimate the V_d of NH₃ in RMNP (Beem et al., 2010; Benedict et al., 2013a, b).

$$\begin{array}{}\text{(2)}& V_{\mathrm{d}}\left({\mathrm{NH}}_{\mathrm{3}}\right)=\mathrm{0.7}\cdot V_{\mathrm{d}}\left({\mathrm{HNO}}_{\mathrm{3}}\right),\end{array}$$

2.5 Bidirectional flux modelling of NH₃

Bidirectional NH₃ fluxes are simulated across the study period using the dry deposition inferential model described in Massad et al. (2010). This model was selected because it estimates both emissions and deposition of NH₃, uses a compensation point framework to capture these complex dynamics, and takes into account rapidly changing micrometeorology. The simulation framework (Fig. 3) accounts for the bidirectional nature of NH₃ fluxes and allows for deposition and emission. The model determines if the flux will be negative (deposition) or positive (emission) based on the relationship between the atmospheric concentration (χ_a) at a given reference height (z) and the compensation point (${\mathit{\chi }}_{z_{\mathrm{0}}}$) at a defined distance (d) above the roughness length (z₀).

Figure 3 Dry deposition inferential model proposed in Massad et al. (2010). The table describes each model element. Arrows next to each flux show the allowed flux directions of the given pathway.

A conceptual diagram of resistances and compensation points is shown in Fig. 3. Aerodynamic (R_a) and laminar boundary layer resistance (R_b) capture the effects of turbulent and diffusive transfer from the atmosphere to the surface, respectively. R_a was calculated according to Thom (1975), where z is 25.35 m, d is 7.15 m, and the roughness length is 1.65 m. The stability functions Ψ_H and Ψ_M for scalars and momentum, respectively, are empirical relationships dependent on L (Thom 1975). Displacement and roughness length were provided from the RMNP NEON Tower (NEON, 2023).

$$\begin{array}{}\text{(3)}& R_{\mathrm{a}}={\left(k\cdot u^{\ast }\right)}^{-\mathrm{1}}\cdot \left(\ln \left({\displaystyle \frac{z-d}{z_{\mathrm{0}}}}\right)-{\mathrm{\Psi }}_{\mathrm{H}}+{\mathrm{\Psi }}_{\mathrm{M}}\right),\end{array}$$

R_b is modeled as described in Xiu and Pleim (2001), where γ_air is the kinematic diffusivity of air, and $D_{{\mathrm{NH}}_{\mathrm{3}}}$ is the diffusivity of NH₃.

$$\begin{array}{}\text{(4)}& R_{\mathrm{b}}={\displaystyle \frac{\mathrm{5}}{u^{\ast }}}\cdot {\left({\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{air}}}{D_{{\mathrm{NH}}_{\mathrm{3}}}}}\right)}^{\mathrm{2}/\mathrm{3}},\end{array}$$

In-canopy resistance (R_g) is the sum of aerodynamic resistance within the canopy (R_ac) and ground boundary layer resistance (R_bg). R_ac was calculated based on Nemitz et al. (2001) using Eq. (5), where α is a height dependent constant calculated using Eqs. (16) and (17) from Massad et al. (2010).

$$\begin{array}{}\text{(5)}& R_{\mathrm{ac}(d+z_{\mathrm{0}})}={\displaystyle \frac{{\mathit{\alpha }}_{(d+z_{\mathrm{0}})}}{u^{\ast }}}\end{array}$$

Ground boundary layer resistance (R_bg) is based on Nemitz et al. (2001), where u_g is the wind speed at the ground, which we approximate as 5 % of the wind speed at tower top (25 m), and z_l is the upper bound height of the logarithmic wind profile above the ground, which we approximate as 10 % of the canopy height (Nemitz et al., 2001).

$$\begin{array}{}\text{(6)}& R_{\mathrm{bg}}=\left({\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{air}}}{D_{{\mathrm{NH}}_{\mathrm{3}}}}}-\ln \left({\displaystyle \frac{D_{{\mathrm{NH}}_{\mathrm{3}}}}{k\cdot u_{\mathrm{g}}\cdot z_{\mathrm{l}}}}\right)\right)\cdot {\displaystyle \frac{\mathrm{1}}{k\cdot u_{\mathrm{g}}}},\end{array}$$

Stomata resistance (R_st) captures the diffusion of NH₃ through plant stomata and is calculated as a minimum value related to the plant type proposed by Hicks et al. (1987). Further parameterization proposed by Nemitz et al. (2001) was used here to calculate R_st, where SR (W m⁻²) is the solar radiation. The minimum value for R_st (225 s m⁻¹) was determined using Table 1 of Zhang et al. (2003), assuming 75 % of the land surface was evergreen needleaf trees and 25 % was deciduous broadleaf trees and shrubs.

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{R}_{\mathrm{st}}& =min\mathit{\left\{}\mathrm{5000}\right(\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}),\mathrm{225}(\mathrm{s}{\mathrm{m}}^{-\mathrm{1}})\\ & \cdot \left(\mathrm{1}+\left({\displaystyle \frac{\mathrm{180}}{SR}}\right)\right)\mathit{\}},\end{array}\end{array}$$

Cuticular resistance (R_w) was calculated according to the proposed corrected parameterization as described in Massad et al. (2010), for a forest ecosystem. When relative humidity (RH) is below 100 %, Eq. (8) is used, and when RH exceeds or is equal to 100 %, Eq. (9) is used.

$$\begin{array}{}\text{(8)}& {\displaystyle }{R}_{\mathrm{w}}=\mathrm{31.5}\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{AR}}}\cdot e^{\left(\mathrm{0.0318}(\mathrm{100}-\mathrm{RH})\right)},\text{(9)}& {\displaystyle }{R}_{\mathrm{w}}={\displaystyle \frac{\mathrm{31.5}}{\mathrm{AR}}},\end{array}$$

In both equations, AR is the acid ratio, which is calculated using the molar ratio of acids and bases in the atmosphere. The calculated acid ratio had a mean value of 1.3, a minimum of 0.22, and a maximum of 11.6. Acid ratios were the largest in the winter months.

$$\begin{array}{}\text{(10)}& \mathrm{AR}={\displaystyle \frac{\mathrm{2}\cdot \left[{\mathrm{SO}}_{\mathrm{2}}\right]+\left[{\mathrm{HNO}}_{\mathrm{3}}\right]}{\left[{\mathrm{NH}}_{\mathrm{3}}\right]}},\end{array}$$

For this study period, the acid ratio was calculated using weekly CASTNET measurements of SO₂ and HNO₃ paired with our measurements of NH₃.

χ_st was calculated according to Massad et al. (2010). In the stomata compensation point (Eq. 11), Γ_st is the emission potential of the stomata and is approximated as 29 based on vegetation samples from the area surrounding the NEON Tower. The sampling methods and determination of this value can be found in the Supplement. Emission potentials describe the potential for surface emission.

$$\begin{array}{}\text{(11)}& {\mathit{\chi }}_{\mathrm{st}}={\displaystyle \frac{\mathrm{2.7457}\times {\mathrm{10}}^{\mathrm{15}}}{T}}\cdot e^{\left(-\frac{\mathrm{10378}}{T}\right)}\cdot {\mathrm{\Gamma }}_{\mathrm{st}},\end{array}$$

χ_g was calculated according to Eq. (3) through Eq. (5) of Stratton et al. (2018). In Eq. (12), TAN is the concentration of total ammoniacal N (the sum of NH₃ and NH${}_{\mathrm{4}}^{+}$) in the soil aqueous phase (mg kg⁻¹), K_H is the Henry constant, and K_a is the equilibrium constant. TAN was estimated at 10.6 mg kg⁻¹ based on soil measurements in RMNP from Stratton et al. (2018). NH₃ flux simulations are very sensitive to the TAN value. The Supplement includes a test of the sensitivity of the flux results to TAN values within one standard deviation for the measurements taken by Stratton et al. (2018).

$$\begin{array}{}\text{(12)}& {\mathit{\chi }}_{\mathrm{g}}={\displaystyle \frac{K_{\mathrm{H}}}{\mathrm{1}+\left({\mathrm{10}}^{-\mathrm{pH}}\right)/\left(K_{\mathrm{a}}\right)}}\cdot \mathrm{TAN},\end{array}$$

K_H and K_a were predicted using Eqs. (13) and (14) based on the models of Montes et al. (2009), where T is temperature.

$$\begin{array}{}\text{(13)}& {\displaystyle }{K}_{\mathrm{H}}=\left({\displaystyle \frac{\mathrm{0.2138}}{T}}\right)\cdot {\mathrm{10}}^{\left(\mathrm{6.123}-\mathrm{1825}/T\right)},\text{(14)}& {\displaystyle }{K}_{\mathrm{a}}={\mathrm{10}}^{\left(\mathrm{0.05}-\frac{\mathrm{2788}}{T}\right)},\end{array}$$

χ_c, Eq. (15) below, was calculated using Eq. (12) from Massad et al. (2010).

$$\begin{array}{}\text{(15)}& {\mathit{\chi }}_{\mathrm{c}}={\displaystyle \frac{\begin{array}{c}{\mathit{\chi }}_{\mathrm{a}}\cdot {\left(R_{\mathrm{a}}\cdot R_{\mathrm{b}}\right)}^{-\mathrm{1}}\\ +{\mathit{\chi }}_{\mathrm{st}}\cdot \left[{\left(R_{\mathrm{a}}\cdot R_{\mathrm{st}}\right)}^{-\mathrm{1}}+{\left(R_{\mathrm{b}}\cdot R_{\mathrm{st}}\right)}^{-\mathrm{1}}+{\left(R_{\mathrm{g}}\cdot R_{\mathrm{st}}\right)}^{-\mathrm{1}}\right]\\ +{\mathit{\chi }}_{\mathrm{g}}\cdot (R_{\mathrm{b}}\cdot R_{\mathrm{g}}{)}^{-\mathrm{1}}\end{array}}{\begin{array}{c}{\left(R_{\mathrm{a}}\cdot R_{\mathrm{b}}\right)}^{-\mathrm{1}}+{\left(R_{\mathrm{a}}\cdot R_{\mathrm{st}}\right)}^{-\mathrm{1}}+{\left(R_{\mathrm{a}}\cdot R_{\mathrm{w}}\right)}^{-\mathrm{1}}\\ +{\left(R_{\mathrm{b}}\cdot R_{\mathrm{g}}\right)}^{-\mathrm{1}}+{\left(R_{\mathrm{b}}\cdot R_{\mathrm{st}}\right)}^{-\mathrm{1}}\\ +{\left(R_{\mathrm{b}}\cdot R_{\mathrm{w}}\right)}^{-\mathrm{1}}+{\left(R_{\mathrm{g}}\cdot R_{\mathrm{st}}\right)}^{-\mathrm{1}}+{\left(R_{\mathrm{g}}\cdot R_{\mathrm{w}}\right)}^{-\mathrm{1}}\end{array}}},\end{array}$$

Compensation point at the displacement height (d) above the roughness length (z₀) is calculated using Eq. (17) below as proposed in Massad et al. (2010). ${\mathit{\chi }}_{z_{\mathrm{0}}}$ takes all other compensation points and resistances into account.

$$\begin{array}{}\text{(16)}& {\mathit{\chi }}_{z_{\mathrm{0}}}={\displaystyle \frac{\left(\frac{{\mathit{\chi }}_{\mathrm{a}}}{R_{\mathrm{a}}}+\frac{{\mathit{\chi }}_{\mathrm{g}}}{R_{\mathrm{g}}}+\frac{{\mathit{\chi }}_{\mathrm{c}}}{R_{\mathrm{b}}}\right)}{\left(\frac{\mathrm{1}}{R_{\mathrm{a}}}+\frac{\mathrm{1}}{R_{\mathrm{g}}}+\frac{\mathrm{1}}{R_{\mathrm{b}}}\right)}},\end{array}$$

Finally, the total flux was calculated following Eq. (17) (Massad et al., 2010). NH₃ flux is calculated in this framework as a difference between the ${\mathit{\chi }}_{z_{\mathrm{0}}}$ and χ_a, scaled by R_a.

$$\begin{array}{}\text{(17)}& F_{{\mathrm{NH}}_{\mathrm{3}}}={\displaystyle \frac{{\mathit{\chi }}_{z_{\mathrm{0}}}-{\mathit{\chi }}_{\mathrm{a}}}{R_{\mathrm{a}}}},\end{array}$$

Total exchange flux ($F_{{\mathrm{NH}}_{\mathrm{3}}}$) from the dry deposition inferential model gives the direction and magnitude of NH₃ fluxes.

3 Results and Discussion

3.1 Simulated bidirectional exchange fluxes of NH₃

Bidirectional fluxes were simulated using the 30 min NH₃ concentration data set and in situ meteorological data as inputs to the Massad et al. (2010) model, described above. NH₃ concentration, ${\mathit{\chi }}_{z_{\mathrm{0}}}$, and fluxes have a strong seasonal cycle in RMNP (see Fig. 4). NH₃ flux direction is determined by the difference between ${\mathit{\chi }}_{z_{\mathrm{0}}}$ and χ_a (Fig. 4a). When NH₃ concentration exceeds the compensation point, NH₃ is deposited to the surface (a negative flux value). Both NH₃ concentrations and deposition fluxes tend to be greatest during the summer (June, July, and August), with 47 % of NH₃ modeled annual dry deposition occurring during June, July, and August. NH₃ fluxes also had the largest variability in the summer. Deposition in the spring (March, April, and May) closely follows, with 43 % of NH₃ modeled annual dry deposition occurring during March, April, and May. During all seasons, there are periods of net emission from the surface (Fig. 4b). The largest periods of net emission occur in the summer. Daily NH₃ emission fluxes are most common in the winter (December, January, and February), although they are typically smaller than deposition fluxes in the spring and summer.

Figure 4 Daily mean values of: (a) Daily mean χ_a and ${\mathit{\chi }}_{z_{\mathrm{0}}}$, and (b) NH₃ flux.

Total modeled NH₃ flux can be broken down into stomatal, ground, and cuticular fluxes. Figure 5 shows the distribution of simulated NH₃ fluxes for each of these components.

Deposition is driven primarily by stomatal and cuticular fluxes, while ground emission fluxes are sometimes observed. Winter periods of net emission (see Fig. 4b) are driven by the ground flux. One potential limitation of the model used for simulations is that it does not consider snow cover on the ground, which could alter winter fluxes in RMNP.

Figure 5 Total NH₃ simulated fluxes are separated into their component fluxes (stomatal, ground, and cuticular). Simulated fluxes are shown for the entire study period. Boxes show the 25th and 75th percentiles, and whiskers are determined at 1.5 times the interquartile range.

NH₃ concentrations at RMNP are impacted by emission and transport patterns, which can both increase daytime NH₃ concentrations. NH₃ emissions from agricultural sources have a strong diel pattern driven by volatilization during warmer daytime temperatures. At RMNP, transport from these regions is driven on many days by the mountain-plains circulation, which begins in the late morning and transports polluted air masses westward and upslope to the park (Gebhart et al., 2011). Previous studies have demonstrated that the upslope transport from sources in the Front Range has impacts on deposition and air concentrations in RMNP (Benedict et al., 2018; Pan et al., 2021). During this study, the largest χ_a values are also observed during upslope transport from source regions in the CO Front Range. These source regions likely disproportionately contribute to NH₃ dry deposition because the difference between χ_a and ${\mathit{\chi }}_{z_{\mathrm{0}}}$ drives the sign and magnitude of the NH₃ flux. On mornings following overnight dew formation, local volatilization from evaporating dew has also been shown to increase morning NH₃ concentrations (Wentworth et al., 2016). This phenomenon was observed in RMNP and corresponds to the increase in the NH₃ diel pattern around 10:00 MT observed in Fig. 6a. One limitation of the bidirectional flux model used is that NH₃ uptake and emission from dew are not simulated. NH₃ concentration, compensation point, and simulated fluxes each have a strong diel pattern, which peaks during the middle of the day (see Fig. 6). The peak value typically occurs close to 13:00. The soil temperature diel pattern contributes to a higher ${\mathit{\chi }}_{z_{\mathrm{0}}}$ during the middle of the day. The annual cycle of soil temperature also contributes to the higher ${\mathit{\chi }}_{z_{\mathrm{0}}}$ observed in summer. Although both NH₃ concentration and compensation point peak during mid-day, we also observe peak deposition fluxes during the middle of the day, indicating that the influence of the diel pattern of NH₃ concentration is stronger than that of the compensation point diel pattern.

Figure 6 Diel pattern of: (a) NH₃ concentration, (b) simulated ${\mathit{\chi }}_{z_{\mathrm{0}}}$, and (c) NH₃ fluxes are shown for the full study period in RMNP. Boxes show the 25th and 75th percentiles, and whiskers are determined at 1.5 times the interquartile range.

To understand the relative importance of NH₃ deposition in RMNP, NH₃ flux simulation results are combined with NADP/NTN wet deposition fluxes and dry deposition fluxes for particulate ammonium (NH${}_{\mathrm{4}}^{+}$) and nitrate (NO${}_{\mathrm{3}}^{-}$) and gaseous HNO₃ derived from CASTNET concentration observations and MLM deposition velocities, to construct an updated seasonal and annual budget of inorganic N deposition at RMNP. This N_r deposition budget for all measured inorganic species is shown in Fig. 7a. Due to the lack of current measurements, wet and dry deposition of organic nitrogen are not included. Benedict et al. (2013b) reported annual organic nitrogen wet deposition of 0.6 kg N ha⁻¹ yr⁻¹ during their 2008–2009 study. NH₃ dry deposition is the net surface flux from the simulations using 30 min NH₃ concentration. The inorganic annual N_r deposition budget totals 3.4 kg N ha⁻¹ yr⁻¹, with the largest contributions coming from NH${}_{\mathrm{4}}^{+}$ wet deposition (1.34 kg N ha⁻¹ yr⁻¹), NH₃ net dry deposition (0.12 kg N ha⁻¹ yr⁻¹), NO${}_{\mathrm{3}}^{-}$ wet deposition (0.71 kg N ha⁻¹ yr⁻¹), and HNO₃ dry deposition (0.33 kg N ha⁻¹ yr⁻¹). Overall, reduced N_r deposition comprises 59 % of the total inorganic N deposition to RMNP. NH₃ dry deposition comprises 4 % of total inorganic N_r deposition. Simulated NH₃ dry deposition (0.11 kg N ha⁻¹ yr⁻¹) is smaller than the value estimated by Benedict et al. (2013b) during their 2008–2009 study (0.66 kg N ha⁻¹ yr⁻¹). The previous value estimated NH₃ dry deposition velocity by scaling the HNO₃ dry deposition velocity by 0.7, instead of simulating the bidirectional exchange of NH₃.

Figure 7 Reactive nitrogen deposition is shown for all species with measured concentrations or deposition for the full year of study. Wet deposition data are from the NADP NTN site at Beaver Meadows. NH₃ dry deposition is modeled using the bidirectional framework from Massad et al. (2010) and 30 min NH₃ concentration data. Dry deposition of HNO₃(g), NH${}_{\mathrm{4}}^{+}$(p), and NO${}_{\mathrm{3}}^{-}$(p) are calculated from the nearby CASTNET site concentration data and deposition velocities from the U.S. EPA MLM. Panel (a) has the annual deposition of all measured species. Panel (b) has deposition of all measured N_r species grouped by month. Only one period of wet deposition was collected by the NADP NTN site during November 2021.

Speciated monthly dry deposition is plotted in Fig. 7b to probe the seasonality of N_r deposition in RMNP. Net dry deposition of NH₃ was largest during May and August. Total inorganic N_r deposition peaked during May, due to increased wet deposition. Reduced N_r deposition exceeded oxidized N_r deposition in October, December, February, March, April, May, July, and August. Excluding November, where only one period of wet deposition was recorded by the NADP NTN site, reduced N_r deposition had a fractional contribution ranging from 43 % to 74 %. In November and January, net NH₃ emission was estimated from the surface.

3.2 Impacts of biweekly NH₃ concentration data on simulated fluxes

The use of low time-resolution NH₃ concentrations for flux simulations can produce a low bias compared with fluxes simulated using higher time-resolution NH₃ concentrations. Simulated NH₃ fluxes have a strong diel pattern when simulated at 30 min resolution (see Fig. 6c), due to changes in NH₃ concentration and meteorology. These complex dynamics are averaged out when an average NH₃ concentration is used, which leads to an underestimation in deposition. Here, we demonstrate that a site-specific correction can be generated to account for the bias introduced by lower time-resolution NH₃ concentration data. Our methods differ from Schrader et al. (2018) in 3 major ways: (i) in situ data is used for both the higher frequency, 30 min NH₃ concentration, and meteorology, (ii) biweekly passive NH₃ data is used instead of monthly NH₃ data, and (iii) Massad et al. (2010) is used as described instead of using a simplified parameterization. The results of the 30 min NH₃ and Biweekly NH₃ bidirectional NH₃ flux simulations are compared to generate a site-specific factor to correct for any low bias in the lower time-resolution flux calculations. Simulated fluxes at biweekly time-resolution (Fig. 8) using the two NH₃ concentration data sets are well correlated (R²=0.88), and the NH₃ flux simulation using biweekly integrated NH₃ data can be corrected to match the control flux simulation using a linear fit (slope = 1.03, y-intercept = −1.689). As noted above, RMNP has few two-week periods of net NH₃ emission, and the efficacy of this method should be confirmed at a location with more extensive periods of net NH₃ emission. In particular, NH₃ fluxes above managed agricultural land could differ significantly from the pattern observed in RMNP. This study also focused on fluxes above a forest canopy, and results could differ for grassland ecosystems, which also occur in RMNP. To determine the efficacy in other locations, future investigations should select several sites with different land surface types and NH₃ concentrations to make biweekly and high time-resolution measurements for a year.

Figure 8 Bidirectional NH₃ flux simulated at 30 min resolution is plotted for 30 min NH₃ concentration data and biweekly integrated NH₃ concentration data. Fluxes are given as net flux over a two-week period. The least squares linear regression is plotted for the data.

Considering the net flux of NH₃ across the full study period, using the best available time-resolution of 30 min, we find a total annual net NH₃ dry deposition flux of 0.11 kg N ha⁻¹ yr⁻¹ (Fig. 9). The estimated NH₃ dry deposition drops by 45 % to 0.06 kg N ha⁻¹ yr⁻¹ using biweekly vs. 30 min NH₃ concentration measurements. The annual NH₃ dry deposition flux increases to 0.78 kg N ha⁻¹ yr⁻¹ when simulating fluxes in a deposition-only unidirectional framework where the NH₃ deposition velocity is scaled as 0.7 times the nitric acid deposition velocity (estimated by the US EPA MLM), an approach previously used for RMNP N deposition budgets (Beem et al., 2010; Benedict et al., 2013a, b).

Figure 9 Annual NH₃ dry deposition at the NEON Flux Tower in RMNP is shown for three bidirectional simulations using three sets of NH₃ concentration data (30 min NH₃, Biweekly NH₃, and Average Diel Pattern NH₃) and one unidirectional simulation. Each simulation was run at 30 min time steps with meteorological parameters from the NEON Flux Tower. The unidirectional NH₃ flux is calculated using biweekly NH₃ concentrations. NH₃ deposition velocities are calculated as 0.7 times the HNO₃ deposition velocity from the U.S. EPA MLM.

Bidirectional flux simulations using biweekly NH₃ data with an average diel pattern of NH₃ yield the same annual NH₃ dry deposition flux as the simulations run using 30 min NH₃ concentration. This indicates that capturing daily variability in NH₃ concentration profiles is not critical to accurately simulating the annual NH₃ flux. Application of an annual averaged diel pattern misses the highest NH₃ concentrations (Fig. 10); however, across a full year of data, the diel pattern effectively captures the net surface flux. Despite the scatter in Fig. 10a, fluxes simulated with an average diel pattern NH₃ data set are well correlated with simulations using 30 min NH₃ concentrations (R²=0.59) and have a fit close to unity. The daily mean fluxes (Fig. 10b and c) of each simulation have similar seasonal patterns, with periods of net emission and deposition aligned between simulations.

Figure 10 NH₃ fluxes simulated with 30 min NH₃ concentrations and annual average diel pattern NH₃ concentrations are shown for the full year of data. Panel (a) directly compares 30 min simulated fluxes for each data set. Panels (b) and (c) show the daily mean fluxes for simulations with 30 min NH₃ concentration and average diel pattern NH₃ concentration, respectively.

At RMNP, there is a large daily variability in concentration due especially to changes in upslope transport. When an air mass arrives from the Colorado Front Range and NE Colorado, NH₃ concentrations rise significantly due to the large emission sources upwind. For the comparison shown in Fig. 10, the diel pattern was determined using a full year of NH₃ concentration data. Fluxes were also simulated using diel patterns determined with only a month of data, to probe the necessary length of measurements to generate an effective diel pattern. Annual deposition from all flux simulations using each different monthly diel pattern fell within 2 % of the annual deposition using the annual average diel pattern. Therefore, in RMNP, one month of 30 min measurements appears sufficient to generate a diel pattern that will effectively correct the net NH₃ surface flux. Other locations may have larger and/or more complex variability in NH₃ diel patterns and may require longer periods of data collection to establish an NH₃ diel pattern.

3.3 Impacts of reanalysis meteorological data on simulated NH₃ fluxes

Bidirectional exchange models require several meteorological and soil parameters, which may not be readily available for many locations of interest. Reanalysis data can provide meteorological inputs for locations where required in situ meteorological and soil measurements are unavailable. To examine the impact on flux simulation accuracy resulting from this substitution at RMNP, the same simulations of NH₃ bidirectional fluxes were run using ERA5 meteorology and soil data. 30 min NH₃ simulations run with reanalysis data inputs are well correlated (R²=0.76) with 30 min NH₃ simulations run with in situ data inputs (see Fig. 11), but overestimate the annual NH₃ deposition flux by a factor of 2. From Fig. 11, we find that the use of ERA5 reanalysis data in the simulation of NH₃ bidirectional fluxes introduces a low bias to the flux magnitude in RMNP compared to in situ meteorological data, for both positive (emission) and negative (deposition) fluxes. However, because the decrease in deposition fluxes is smaller than the decrease in emission fluxes, we observe an annual overestimation from simulations using ERA5.

Figure 11 Bidirectional NH₃ flux simulated with ERA5 meteorology and NEON meteorology at 30 min resolution using the 30 min NH₃ concentration. The least squares linear regression is plotted for the data in red.

The low bias for fluxes simulated using ERA5 reanalysis data is investigated further to explore what parameters influence this bias. Net NH₃ fluxes are simulated using Eq. (17), which relies on ${\mathit{\chi }}_{z_{\mathrm{0}}}$, NH₃ concentration, and R_a. We find that the simulations using reanalysis data generate ${\mathit{\chi }}_{z_{\mathrm{0}}}$, which agree well with the simulations using in situ measurements (slope = 1.03, R²=0.98).

Figure 12 Aerodynamic resistances are shown for simulations using in situ meteorological data from the NEON flux tower and reanalysis meteorological data from ERA5. The diel patterns are shown in (a) and (b), respectively. Panel (c) directly compares simulated R_a values using NEON in situ and ERA5 reanalysis data.

Although the general diel pattern of R_a is well captured using reanalysis data, R_a magnitudes differ substantially between the two simulations (Fig. 12a and b), with the largest difference occurring overnight. Maximum R_a values from the reanalysis simulations are an order of magnitude larger than those derived using in situ meteorology. A comparison of the two data sets shows (Fig. 12c) a typical enhancement of approximately a factor of four. Increased R_a values result in lower simulated NH₃ fluxes. The R_a bias is likely driven by differences in the u_∗ and L, which are used to calculate R_a. ERA5 data underestimates u_∗ by a factor of 5 when compared with the in situ NEON data (slope = 0.2). The in situ NEON data also sets a minimum u_∗ value (0.2 m s⁻¹), while the ERA5 data allows u_∗ values below 0.2 m s⁻¹. Comparisons of all meteorological parameters used can be found in the Supplement. This discrepancy in modeled R_a may be due to the gridded nature of reanalysis data, which represents a large area of variable land types and complex topography using only a single value (Hogrefe et al., 2023). Previous work has identified heat and moisture fluxes as large areas of uncertainty in ERA5 Reanalysis (Kong et al., 2022; Mayer et al., 2022). Two case studies were conducted to probe the relative importance of u_∗ and L. The case studies are described in the Supplement. Differences in R_a were impacted by both u_∗ and L, accounting for 23 % and 10 %, respectively, of the discrepancy between in situ and ERA5 flux simulations.

4 Conclusion

Fluxes of NH₃(g) can be simulated using a bidirectional model, which uses rapidly changing meteorology paired with air concentrations and soil parameters to infer flux direction and magnitude. We use a bidirectional NH₃ flux model, proposed by Massad et al. (2010), to simulate a year of NH₃ fluxes above a subalpine forest ecosystem in Rocky Mountain National Park. The net NH₃ dry deposition to the ecosystem is estimated at 0.11 kg N ha⁻¹ yr⁻¹, comprising 4 % of total inorganic reactive nitrogen deposition. This is significantly lower than previous estimates for RMNP, which did not consider the bidirectional nature of the exchange. Due to the observed low bias in passive NH₃ observations and the sensitivity of simulations to NH₃ concentrations, this is likely a low bound. The sensitivity of NH₃ flux modelling to χ_a was tested by scaling the input concentration by 9 % to account for the error discussed in Puchalski et al. (2011). This resulted in an annual deposition increase of 47 %, indicating the importance of accurate NH₃ measurements for flux modelling. Additionally, since the highest NH₃ concentrations occur during upslope events, the sources contributing to these events likely have a disproportionate effect on deposition. One limitation of this model is the exclusion of snow cover, which could significantly change NH₃ fluxes in the winter, when RMNP has frequent snow events. To probe the impact of snow cover, a sensitivity test was conducted setting χ_g equal to zero during the winter (December, January, and February), which increased annual deposition by 0.06 kg N ha⁻¹ yr⁻¹. However, this analysis does not take into account how the surface differences may change NH₃ fluxes above snow. Future work should investigate NH₃ fluxes above snow cover to better simulate the exchange of NH₃ in regions with snow.

Due to the cost and technical challenges of making continual, high time-resolution NH₃ concentration measurements, there is growing interest in using integrated biweekly passive NH₃ measurements, such as those from the NADP AMoN network, for flux simulations. Here, we establish that a site-specific correction can be used to correct a bias introduced by using lower time-resolution passive NH₃ measurements over the studied forest canopy in RMNP. We also establish that an average NH₃ diel pattern can be used to interpolate 30 min NH₃ concentration and correct for the bias introduced by passive NH₃ measurements. In RMNP, a month of measurements proved sufficient to determine the diel pattern used for flux simulations. The correction factor and diel pattern, however, likely vary by location due to differences in ecosystem characteristics and factors influencing NH₃ concentrations. Due to the potential regional differences and changes associated with land surface type, additional sites should be studied to assess the impact of measurement time-resolution on NH₃ flux simulations. To understand the seasonal variability in diel pattern and efficacy of diel pattern application for flux simulations, measurements should be conducted for a full year.

Local micrometeorological and soil measurements are also frequently unavailable, making the use of reanalysis data a desirable alternative for NH₃ flux simulations. In our location, the use of reanalysis data adds a bias that leads to overestimates of net NH₃ deposition. We found it was possible to apply a correction to address this bias, but this factor likely varies by location, in particular over different land surface types within a reanalysis grid cell. Future studies should explore the relationship between in situ measurements and reanalysis products above different land surface types, varied topography, and in different regions. Understanding how to correct for the biases introduced through the use of reanalysis data would allow improved modelling of NH₃ bidirectional fluxes in regions lacking high time-resolution measurements.

In this analysis, we simulated the bidirectional exchange of NH₃ above a forest ecosystem using the model proposed in Massad et al. (2010). However, there are other bidirectional exchange models (e.g., Zhang et al., 2010; Pleim et al., 2013) and their simulated fluxes may differ significantly from the model used here (Jongenelen et al., 2025). In the bidirectional exchange model used here, we observe that the selected inputs for NH₃ concentration and meteorological data may introduce biases into the simulated NH₃ fluxes. This may also be true for the other models when simulating NH₃ bidirectional exchange, a good topic for future research.