Source regions contributing to excess reactive nitrogen deposition in the Greater Yellowstone Area (GYA) of the United States

Research has shown that excess reactive nitrogen (N_r) deposition in the Greater Yellowstone Area (GYA) of the United States has passed critical load (CL) thresholds and is adversely affecting sensitive ecosystems in this area. To better understand the sources causing excess N_r deposition, the Comprehensive Air Quality Model with Extensions (CAMx), using Western Air Quality Study (WAQS) emission and meteorology inputs, was used to simulate N_r deposition in the GYA. CAMx's Particulate Source Apportionment Technology (PSAT) was employed to estimate contributions from agriculture (AG), oil and gas (OG), fire (Fire), and other (Other) source sectors from 27 regions, including the model boundary conditions (BCs) to the simulated N_r for 2011. The BCs were outside the conterminous United States and thought to represent international anthropogenic and natural contributions. Emissions from the AG and Other source sectors are predominantly from reduced N and oxidized N compounds, respectively. The model evaluation revealed a systematic underestimation in ammonia (NH₃) concentrations by 65 % and overestimation in nitric acid concentrations by 108 %. The measured inorganic N wet deposition at National Trends Network sites in the GYA was overestimated by 31 %–49 %, due at least partially to an overestimation of precipitation. These uncertainties appear to result in an overestimation of distant source regions including California and BCs and an underestimation of closer agricultural source regions including the Snake River valley. Due to these large uncertainties, the relative contributions from the modeled sources and their general patterns are the most reliable results. Source apportionment results showed that the AG sector was the single largest contributor to the GYA total N_r deposition, contributing 34 % on an annual basis. A total of 74 % of the AG contributions originated from the Idaho Snake River valley, with Wyoming, California, and northern Utah contributing another 7 %, 5 %, and 4 %, respectively. Contributions from the OG sector were small at about 1 % over the GYA, except in the southern Wind River Mountain Range during winter where they accounted for more than 10 %, with 46 % of these contributions coming from OG activities in Wyoming. Wild and prescribed fires contributed 18 % of the total N_r deposition, with fires within the GYA having the highest impact. The Other source category was the largest winter contributor (44 %) with high contributions from California, Wyoming, and northern Utah.

1 Introduction

The Greater Yellowstone Area (GYA) (see Fig. 1) of the United States, with Yellowstone National Park (YNP) and Grand Teton National Park (GTNP) at its core, is one of the largest remaining intact ecosystems in the northern temperate zone and features diverse wildlife, alpine lakes, forests, and geologic wonders (Keiter and Boyce, 1994; NPS, 2017). Increasing concentrations of reactive nitrogen (N_r) compounds in air, rain, and snowpack samples over the GYA have been reported in the past 30 years and linked to N_r emissions from human activities (Clow et al., 2003; Blett et al., 2011; IMPROVE, 2011; Sullivan et al., 2011; USGS, 2014; NADP, 2016; Nanus et al., 2017; also, see Fig. S1). The inorganic wet N_r deposition rates measured at high-elevation National Trends Network (NTN) sites within the GYA in 2010 were 2.5–3.5 kg N ha⁻¹ yr⁻¹, compared with 1.5–2.5 kg N ha⁻¹ yr⁻¹ in 2000 (NADP, 2016). This is relevant to the long-term conservation of the area because as N_r deposition levels increase, they can cross critical load (CL) thresholds, at which negative effects to sensitive ecosystem components can occur (Porter et al., 2005; Pardo et al., 2011). Additional concerns posed by enhanced N_r deposition include lake acidification, loss of lichen biodiversity, and eutrophication (Baron, 2006; Blett et al., 2011; NADP, 2016). While ecosystem changes due to excess N_r deposition over Class I areas including the GYA have been documented (e.g., Baron et al., 2011; Saros et al., 2011; Sullivan et al., 2011; Spaulding et al., 2015; Nanus et al., 2017), the origins, chemical composition, and spatial and temporal changes in the deposition over this region are not as well understood.

Total N_r is a mix of oxidized and reduced inorganic nitrogen (N) and organic N compounds that are chemically and biologically active in the Earth's biosphere and atmosphere and deposited through wet and dry processes. These compounds arise from a variety of sources, with inorganic oxidized N primarily emitted as nitrogen oxides (NO_x) from fossil fuel combustion, with approximately 25 % from power plants, 50 % from automobiles, and 10 % from other mobile sources, based on annual county-level estimations (EPA, 2015a). Atmospheric reactions of NO_x result in nitric acid (HNO₃), particulate nitrate (PNO₃), and other compounds. Reduced N arises primarily from ammonia (NH₃) gas emissions from agricultural activities, which can react with acidic aerosols to form ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) compounds (Galloway et al., 2004). Mobile sources are also an important source of NH₃ and can be the primary emitter in urban areas (Sun et al., 2014, 2017). Emissions from this sector have large uncertainties and a recent study suggests that on-road NH₃ emissions in the 2011 National Emissions Inventory (NEI) were underestimated by a factor of 2.9 (Fenn et al., 2018). There are hundreds of organic N compounds, including reduced (e.g., amines) and oxidized forms (e.g., alkyl nitrates). Sources of organic N are less well known, but increasing evidence suggests that biomass burning and agriculture (AG) are significant contributors, as are atmospheric reactions of NO_x with volatile organic compounds (Cape et al., 2011; Reay et al., 2012). With the steady decline of NO_x emissions in the United States during past decades as a result of the implementation of the Clean Air Act, the importance of reduced N to the total N deposition budget has increased (Li et al., 2016). Specific to the GYA, local anthropogenic emissions are small, but upwind sources, including agricultural activities in the Snake River valley and northern Utah, wildfires throughout the western United States, energy development in the Upper Green River basin, and anthropogenic activities at urban centers such as Salt Lake City, are larger and likely to be significant contributors to regional N emissions (Prenni et al., 2014).

To better understand the levels and composition of the N_r compounds deposited in the GYA and to help guide strategies to reduce N deposition, the National Park Service (NPS) initiated the Grand Teton Reactive Nitrogen Deposition Study (GrandTReNDS), which included spatially and temporally detailed measurements of N compounds during April to September 2011 (Benedict et al., 2013a; Prenni et al., 2014). It was found that during summer months at the high-elevation sites (e.g., Grand Targhee; see Fig. 2), 62 % of the N deposition was from reduced N and about equally split between dry and wet deposition, and oxidized N only accounted for 27 % of the N deposition budget, with the remaining N in the form of wet-deposited, organic N. Study findings indicate a significant west-to-east gradient in atmospheric NH₃ concentrations, with higher concentrations west of the Teton mountain range. Concurrently measured concentrations of HNO₃, PM_2.5 (particulate matter with an aerodynamic diameter of less than 2.5 µm) nitrate, and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ showed relatively small west-to-east gradients inside GTNP (Benedict et al., 2013a; Prenni et al., 2014).

The origins of N_r transported to the GYA and other remote locations in the western United States have been examined in past modeling studies. Back trajectory analyses have shown that air mass transport to GTNP is predominantly from the west through the Snake River valley and from the southwest through northern Utah (Prenni et al., 2014). Zhang et al. (2012) applied the global chemical transport model (CTM) GEOS-Chem (Bey et al., 2001) using zero-out sensitivity simulations and found that in 2006 natural sources, including lightning and wildfires, contributed more than 10 % of the total N_r deposition over the Teton area. Lee et al. (2016) used the adjoint version of GEOS-Chem to quantify the sources of N_r deposition in eight selected federal Class I areas in 2010 and found a nonnegligible footprint (>20 %) of N_r deposition in the western United States, including GTNP and Rocky Mountain National Park (RMNP), attributed to long-range transport from sources in California, especially during summertime. Mobile NO_x and livestock NH₃ were also found to be major sources of N_r deposition in GTNP. Similar modeling studies focusing on RMNP also suggested the important contributions of distant sources including those from California and other counties and the fact that the contributions from sources of reduced N_r were larger than those from sources of oxidized N_r (Thompson et al., 2015; Malm et al., 2016).

In this work, we add to the growing body of N_r modeling source apportionment studies by conducting a detailed analysis using the Particulate Source Apportionment Technology (PSAT) module within the CAMx (Comprehensive Air Quality Model with extensions) (ENVIRON, 2014) CTM to quantify the seasonal contributions from different source regions and source sectors to N_r throughout the GYA. Compared with previous N_r deposition simulation studies in United States, this work uses tagged reactive tracers to attribute the contributions from four designated emission sectors and 27 designated emission regions to N_r deposition in the GYA with a much higher horizontal grid resolution (12 km) and an up-to-date emission inventory instead of using a zero-out approach (e.g., Zhang et al., 2012) or an adjoint model (e.g., Lee et al., 2016). The model simulation of N_r and its constituents were first evaluated against routine measured data as well as the unique data measured during the GrandTReNDS campaign period (Benedict et al., 2013a; Prenni et al., 2014). N_r deposition from CAMx simulations was also compared with total deposition maps (TDEPs), which were developed for deposition trend analysis and ecological impact assessment (Schwede and Lear, 2014). The detailed source apportionment results are presented here, focusing on seasonal variations and the relative importance to CL exceedance in sensitive ecosystems within the GYA. The discussion of identified model bias and uncertainties in the interpretation of source apportionment results, including the model lateral boundary conditions, the impact of model precipitation to wet deposition simulation, and the impact of ammonium dry deposition velocity on dry deposition are also presented.

Figure 1 Source region partition for CAMx PSAT simulation in this study. The 27 tagged regions are (1) NW Colorado; (2) NE Colorado; (3) SE Colorado; (4) SW Colorado; (5) Upper Green River, Wyoming; (6) Jackson, Wyoming; (7) eastern Wyoming; (8) western Wyoming; (9) Yellowstone; (10) northern Idaho; (11) Snake River valley, Idaho; (12) northern Utah; (13) southern Utah; (14) Nevada; (15) Montana; (16) Washington; (17) Oregon; (18) California; (19) Mexico; (20) New Mexico; (21) Arizona; (22) Texas & Oklahoma; (23) Canada; (24) North Dakota; (25) the Pacific; (26) the far east US; and (27) South Dakota, Kansas, and Nebraska.

Figure 2 Annual NH₃ emission for the 12 km inner modeling domain at focused tagged regions (see Table S2 and Fig. 1 for the details of the 27 source region partition) as well as locations of the monitoring sites at different networks (^a Ammonia Monitoring Network; ^b Clean Air Status and Trends Network; ^c Grand Teton Reactive Nitrogen Deposition Study; ^d Interagency Monitoring of Protected Visual Environments; ^e National Trends Network) used in the model performance evaluation of CAMx nitrogen species concentration and dry–wet deposition in the GYA (the black boundary line). The numbers in the figure are locations for the three sampling sites during GrandTReNDS and the eight Class I areas within the area: (1) Driggs, (2) Grand Targhee, (3) NOAA climate station, (4) Grand Teton National Park, (5) John D. Rockfeller Jr. Memorial Parkway, (6) Yellowstone National Park, (7) Teton Wilderness, (8) Washakie Wilderness, (9) North Absaroka Wilderness, (10) Fitzpatrick Wilderness, and (11) Bridger Wilderness.

2 Modeling system for N_r source apportionment

Modeling simulations for 2011 were conducted using the CAMx version 6.10 (ENVIRON, 2014) with two nested grids. The outer domain (36 km) covered the contiguous United States (CONUS), as well as portions of Canada and Mexico, while the inner domain (12 km) encompassed the western United States and focused on states within the Western Regional Air Partnership (WRAP) (see Fig. 1).

The hourly meteorological inputs for 2011 were generated by the Weather Research and Forecasting (WRF) model (WRF-ARW, version 3.5.1) (Skamarock et al., 2008) and were obtained from the Intermountain West Data Warehouse (IWDW) (http://views.cira.colostate.edu/tsdw/, last access: 30 August 2018). This meteorological simulation performed in comparison to other recent prognostic model applications used in air quality planning (UNC-Chapel Hill and ENVIRON, 2014a).

The emission inventory used by CAMx was primarily derived from the 2011 NEI version 2 (NEI2011v2) (EPA, 2015b) with the Sparse Matrix Operator Kernel Emissions (SMOKE) processing system version 3.0 (Houyoux et al., 2002) for anthropogenic emissions, the Model of Emissions of Gases and Aerosols from Nature (MEGAN) version 2.10 (Guenther et al., 2012) for biogenic emissions, and the WRAP Windblown Dust Model (WRAP-WBD) to estimate wind-driven dust emissions (UNC-Chapel Hill and ENVIRON, 2014b). Emissions from the oil and gas (OG) sector were further updated by the IWDW to represent the best-available inventory for OG activity in the western United States at the time of modeling (UNC-Chapel Hill and ENVIRON, 2014b). The emissions for fire activities (Fire) include agricultural fires, prescribed fires, and wildfires and were generated by the Particulate Matter Deterministic and Empirical Tagging and Assessment of Impacts on Levels (PMDETAIL) study (Moore et al., 2012). PMDETAIL developed 2011 fire emissions using satellite data, ground detection, and burn scar and estimated the plume rise, depending on fire size and type. The hourly, nonsurface fire emissions were allocated to the proper CAMx vertical layers based on the model-predicted planetary boundary layer (PBL) height and the spanning of the plume top and bottom above the ground (Mavko and Morris, 2013).

The boundary conditions for the 36 km domain were estimated from a 2011 global model run using the Model for Ozone and Related chemical Tracers (MOZART) version 4.6 (Emmons et al., 2010). The simulation year of 2011 was preceded by 15 days of spin-up time to minimize the effects of initial conditions. A more-detailed description of the WRF–SMOKE–CAMx modeling platform applied in this study is summarized in Table S1 in the Supplement as well as the 2011 Three-State Air Quality Study (3SAQS) (UNC-Chapel Hill and ENVIRON, 2014b).

For the source apportionment estimates, 27 source regions (Fig. 1), as well as the lateral boundary conditions (BCs), were “tagged” in the CAMx PSAT simulation. In addition, the emissions for each region were further subdivided into four source sectors: (1) AG, (2) OG, (3) Fire, including wildfires and prescribed fires, and (4) the remaining sources labeled as Other. The Other source sector primarily comes from mobile and large point sources, with smaller contributions from natural sources such as lightning. Table S2 provides the annual NH₃ and NO_x emissions used in this modeling study with a breakdown by tagged source regions and source sectors. Figure 2 provides the annual emissions of NH₃ in the inner 12 km domain as well as the monitoring sites or receptor areas used for the model evaluation and analysis. For NH₃ emissions, the AG sector contributed 84.1 % of the total emissions within the 12 km domain, while the OG, Fire, and Other sectors contributed 0.1 %, 4.5 %, and 11.4 %, respectively (Table S2). In the Snake River valley, the AG sector emissions dominate the emission budget. For NO_x emissions, the contribution rankings from the four tagged emission sources are Other (83.8 %), OG (12.8 %), Fire (3.2 %), and AG (0 %). The regions were selected to highlight important source sector contributions to N_r deposition in the GYA. For example, the state of Wyoming was partitioned into five regions (YNP, Jackson, Upper Green River, eastern Wyoming, and western Wyoming) to differentiate the possible source impacts from urban activity in Jackson from energy development in southwestern Wyoming (Blett et al., 2011; NPS, 2017). Significant agricultural operations in the Snake River valley in Idaho, northern Utah, and northeastern Colorado were tagged due to their high ammonia emissions (see Fig. 2) associated with fertilizer application and confined animal feeding operations (Fenn et al., 2003; Clarisse et al., 2009; Prenni et al., 2014). Lastly, wildfires are episodic events (http://wrapfets.org/map.cfm, last access: 30 August 2018) that can have large intermittent contributions to N_r deposition, but they can mask important contributions from other sources that are significant in nonfire years.

Table 1 CAMx model performance for nitrogen species concentrations as well as nitrogen dry–wet depositions evaluated at sites in the AMoN, CASTNet, IMPROVE, and NTN networks as well as the three sites during the GrandTReNDS campaign over the GYA region (see Fig. 1 for site locations) in 2011.

¹ AMoN samples are collected for 2 weeks. ² CASTNet samples are collected for 1 week. ³ GrandTReNDS samples are collected daily. ⁴ IMPROVE 24 h samples are collected every 3 days. ⁵ NH_x= NH₃+PNH₄. ⁶ NTN samples are collected for 1 week. ⁷ Average observation. ⁸ Average simulation. ⁹ Number of sites. ¹⁰ Number of samples; the values in the parentheses are the percentage of valid samples used for model performance evaluation. ¹¹ Pearson's correlation coefficient. ¹² Normalized mean bias. ¹³ Normalized mean error. ¹⁴ Fractional bias. ¹⁵ fractional errors.

CAMx–PSAT treats nitrogen-containing compounds as one of seven species: gaseous NH₃; particulate ammonium (PNH₄); reactive gaseous nitrogen (RGN), which includes primary emissions of NO_x, nitrous acid (HONO), nitrate radical (NO₃), and dinitrogen pentoxide (N₂O₅); gaseous nitric acid (HNO₃); gaseous peroxy nitrogen (TPN), including peroxyacetyl nitrate (PAN) and peroxynitric acid (PNA); gas-phase organic nitrate (NTR); and particulate nitrate (PNO₃). PSAT maintains the source-group identity (i.e., source region and source sector) by apportioning the secondary species to the precursor emissions (ENVIRON, 2014). In the source apportionment comparison results, we report the reduced N_r deposition as the sum of NH₃ and PNH₄ and the oxidized N_r deposition as the sum of RGN, HNO₃, PNO₃, TPN, and NTR in units of kg N ha⁻¹.

3 Evaluation of CAMx-simulated N_r concentration and deposition rates

Acceptable model performance of the regional air quality modeling system is a prerequisite for a credible source apportionment interpretation (Boylan and Russell, 2006; EPA, 2014; Emery et al., 2017). In this work, the CAMx simulation was extensively evaluated against routine monitoring data as well as data collected in the GrandTReNDS special field study (Benedict et al., 2013a; Prenni et al., 2014) and against the nitrogen deposition estimates from the National Atmospheric Deposition Program (NADP, http://nadp.slh.wisc.edu/, last access: 30 August 2018) TDEP hybrid modeling results (Schwede and Lear, 2014). Performance metrics recommended by the EPA's modeling guidance for ozone, PM_2.5, and regional haze attainment demonstrations (Yu et al., 2006; EPA, 2014) were used (see Table 1).

The variables and routine monitoring networks used in the model evaluation were NH₃ concentrations from the Ammonia Monitoring Network (AMoN) (http://nadp.sws.uiuc.edu/AMoN/, last access: 30 August 2018); nitric acid (HNO₃), PNO₃, and PNH₄ concentrations as well as estimated dry deposition fluxes from the Clean Air Status and Trends Network (CASTNet) (https://www.epa.gov/castnet, last access: 30 August 2018); PNO₃ and PNH₄ concentrations from the Chemical Speciation Network (CSN) (https://www3.epa.gov/ttnamti1/speciepg.html, last access: 30 August 2018); PNO₃ concentrations from the Interagency Monitoring of Protected Visual Environments (IMPROVE) network; and wet-deposited inorganic oxidized (${\mathrm{NO}}_{\mathrm{3}}^{-}$) and reduced (${\mathrm{NH}}_{\mathrm{4}}^{+}$) nitrogen and associated precipitation rates from the NADP NTN sites. Each network had a unique sampling frequency and duration (Table 1). The hourly CAMx outputs were aggregated to match the timescales of the measured data. All measurement data flagged as questionable, either due to mal-operation or to insufficient samples to calculate representative values, were excluded from the analysis. Table 1 reports the percentage of valid measurements used for statistical analysis during evaluation time. For most of the N species, the percentage of valid samples is more than 80 %. In general, the CASTNet, IMPROVE, AMoN, and NADP networks sample in rural areas, while the data from the CSN network primarily represent the air quality in urban and suburban settings. Although organic N species were also measured in the GrandTReNDS campaign, we focus on the inorganic N budget comparison, given the large uncertainties for organic N prediction (Jickells et al., 2013) and its incomplete treatment in the model's chemical mechanism. For example, the modeling system does not account for primary emissions of organic N compounds but does include the formation of organic N from the alkylperoxy radical and secondary alkoxy radical (ENVIRON, 2014).

3.1 Evaluation against data in the GYA

The 3SAQS study performed photochemical grid modeling using the same modeling platform and input files as this study (UNC-Chapel Hill and ENVIRON, 2014b) and evaluated the model performance for the western United States. A subset of these results is presented in the Supplement for reference. Model performance statistics for the N species within the GYA area at AMoN, CASTNet, IMPROVE, and NTN network sites (Fig. 1) at different periods in 2011 are presented in Table 1. The biases at the GYA sites are similar to those throughout the west (Table S1) in that the CAMx simulation significantly overestimated HNO₃ with a normalized mean bias (NMB) of 108 % and significantly underestimated NH₃ concentrations with a NMB $=-\mathrm{65}$ %. While the model had skill in reproducing the daily variation in HNO₃, with a correlation coefficient of r=0.71, it had little skill for NH₃, with r=0.2. The overestimation of HNO₃ has also been reported in other regional-scale modeling simulations over the United States (e.g., Baker and Scheff, 2007; Foley et al., 2010; Thompson et al., 2015) with the carbon bond mechanism used in this study. The possible reason for the overestimation of HNO₃ may be due to the uncertainty for the N₂O₅ uptake coefficient setting for heterogeneous reactions (Foley et al., 2010). The poor NH₃ results may be related to the high uncertainty in the NH₃ emission inventory (Clarisse et al., 2009) and important missing physical mechanisms in the model, including the lack of bidirectional NH₃ deposition (Zhang et al., 2010; Bash et al., 2013; Zhu et al., 2015). The GYA area is located downwind of the major agriculture sources in the Snake River valley and northern Utah (Table S2). The incorporation of the bidirectional NH₃ flux mechanism in the model should increase ambient NH₃ concentrations in the GYA and thus decrease the large model underestimation of NH₃ concentrations.

For PNO₃ and PNH₄ simulations in the GYA, CAMx overestimated both species, with better performance for PNH₄ than PNO₃ (3 % vs. 37 %, respectively, in terms of NMB) and better agreement for PNO₃ at CASTNet sites vs. IMPROVE sites (37 % vs. 58 % for NMB, respectively). The errors and biases in the dry deposition fluxes compared to CASTNet values follow the same patterns as in the ambient concentrations, but it should be noted that CASTNet and CAMx use different algorithms to estimate dry deposition velocities, and these model-to-model discrepancies will manifest themselves in the performance evaluations.

Wet deposition measurements from the five NTN sites with sufficient data were available from within the GYA (Fig. S2). Comparisons to CAMx showed that the model captured the general trends in these data with r∼0.32–0.34 but were somewhat biased, with a NMB = 31 % for ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and NMB = 49 % for ${\mathrm{NH}}_{\mathrm{4}}^{+}$. The precipitation simulations were consistently 100 %–200 % higher than the rain gauge measurements at the NTN sites, showing that WRF overestimated the frequency and intensity of precipitation events over the GYA in 2011 (Table 1). However, note that 2011 was a large snowpack year; by May, much of the GYA had 100 %–180 % of the normal snow weather equivalent (USGS, 2014). Precipitation measurements tend to be low during high-snow events.

Figure 3 Model performance for (a, b) seasonal average N_r concentration, (c, d) seasonal accumulated N_r deposition budget, and (e, f) seasonal accumulated precipitation amount at collocated location sites (YNP and Pinedale) over the GYA in 2011. ¹ Clean Air Status and Trends Network. ² Comprehensive Air Quality Model with extensions. ³ National Trends Network. ⁴ Parameter-elevation Regressions on Independent Slopes Model. ⁵ Weather Research and Forecasting model.

Figure 4 Inorganic nitrogen deposition budgets as an absolute (a) and as a percentage (c) as well as precipitation (e), measured at the three core sites during the GrandTReNDS study period (April to September in 2011) with corresponding CAMx simulations (b, d, e). ¹ Grand Teton Reactive Nitrogen Deposition Study. ² Comprehensive Air Quality Model with extensions. ³ Weather Research and Forecasting model.

The seasonal, simulated ambient concentrations and deposition rates are compared against measured CASTNet and NADP data at the YNP and Pinedale monitoring sites in Fig. 3. Seasons refer to winter (December, January, February, DJF), spring (MAM), summer (JJA), and fall (SON). The significant overestimation of HNO₃ is evident in all seasons. Also evident is the poor simulation of the seasonality in N_r deposition, primarily due to the poor reproduction of wet deposition, which is at least partly due to the large errors in the simulated precipitation.

Table S3 provides a comparison of regional CTM performance evaluations against measured N-containing species over the United States from peer-reviewed studies in recent years (e.g., Simon et al., 2012; Bash et al., 2013; Zhang et al., 2013; Yu et al., 2014; Thompson et al., 2015; Li et al., 2017). The model performance results in this study are comparable to these past studies, including the overestimation of HNO₃ and underestimation of NH₃. Resolution of these biases requires additional research and these biases need to be taken into account when interpreting the source attribution of N_r deposition within the GYA.

Figure 5 Annual nitrogen deposition budgets as an absolute (a) and as a percentage (b) as well as annual precipitation amounts (c) from the NADP Total Deposition Map (TDEP) and corresponding CAMx (Comprehensive Air Quality Model with extensions) and WRF (Weather Research and Forecasting model) simulation results in 2011 at eight Class I areas across the GYA (the receptor sites on the x axis are arranged from west to east in the GYA; see Fig. 2). The reported CAMx dry and wet N_r deposition values at the eight Class I areas are the average of the simulation values at corresponding grid cells for each area.

3.2 Evaluation against GrandTReNDS data

The GrandTReNDS campaign provides a unique opportunity to evaluate the capability of CAMx to simulate the N_r compounds and deposition budget. Detailed measurements, including NH₃, were made at three sites that crossed GTNP from west to east: Driggs, in the foothills just west of GTNP (43.74^∘ N, −111.87^∘ W, elevation 1947 m); Grand Targhee, an upper-elevation site on the western edge of GTNP (43.78^∘ N, −110.94^∘ W, elevation 2722 m); and the National Oceanic and Atmospheric Administration (NOAA) climate station site on the eastern edge of GTNP (43.66^∘ N, −110.71^∘ W, elevation 1978 m) (also see Fig. 2). Figure 4 presents the monthly deposition budgets for these three sites during the sampling periods, and Table 1 provides the model performance statistics for the N species concentration and deposition. As shown, the simulation does a poor job of reproducing the total N_r deposition rates both in the month-to-month variation and across the sites. The difference in the dry NH₃ deposition monthly variation between measurements and simulation is mainly due to the difference in associated dry deposition velocity used for calculation. However, consistent with the observations, the simulation shows that wet deposition is larger than dry deposition and that the contribution from reduced N deposition was larger than from the oxidized N deposition at all three sites, although the observed range of 70 %–80 % reduced N was more than the 55 %–68 % simulated in CAMx. The primary cause of this bias was the overestimation in the HNO₃ dry deposition rates, which were 2–3 times larger than those derived from the measured data. This is consistent with the systematic overestimation of HNO₃ concentrations (NMB = 106 % in Table 1). Other biases also exist, including an underestimation in the NH₃ dry deposition, which was somewhat balanced by an overestimation in the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ wet deposition (NMB = 60 %). The underestimation of NH₃ concentration still existed (NMB = −16 %), and one of the possible reasons may be due to the overestimation of HNO₃ in the model pushing excessive partitioning of NH₃ into the particle phase, which can be shown by the better model performance for NH_x (NH_x = NH₃ + PNH₄) simulation (NMB $=-\mathrm{7}$ %) without splitting the gas-particle partition bias.

An additional challenge that affected model performance was the difficulty in estimating precipitation rates. This is shown in Fig. 4, in which the simulated precipitation rates do not reproduce the month-to-month variation and generally were highly overestimated. For example, on average the simulated precipitation at Driggs was more than double the measured precipitation, and it was more than a factor of 4 higher at the NOAA climate station site.

3.3 Evaluation against NADP TDEP

TDEP maps (Schwede and Lear, 2014) are widely used in the land management community to assess total N_r deposition throughout the United States and estimate the critical load exceedances in sensitive ecosystems (Saros et al., 2011; Nanus et al., 2017). TDEP employs a hybrid approach to integrate measurements from multiple networks, including CASTNet and NTN, with Community Multiscale Air Quality (CMAQ) modeling (Byun and Schere, 2006) results for deposition velocities and unmeasured species' dry deposition, as well as PRISM (Parameter-elevation Regressions on Independent Slopes Model) (Daly et al., 1994) high-resolution precipitation estimates for mapping total deposition in the United States (Schwede and Lear, 2014). Both the CAMx simulation in this study and the TDEP results are derived from model simulations and subject to similar errors in emissions and physical and chemical processes. However, with the incorporation of measured wet N_r deposition and N concentration data into the TDEP results, they are expected to be less biased than the deposition results from a purely CAMx simulation.

The TDEP total N_r deposition and the CAMx 2011 simulation in this work exhibited similar spatial and temporal patterns across the western United States; for example, both sets of results show high N_r deposition in the Snake River valley, northern Utah, and across the Wyoming state border area near GTNP, with values >5 kg N ha⁻¹ yr⁻¹. Within the GYA (Fig. S3), the CAMx simulation had higher dry N_r deposition, which was more spatially heterogeneous than the corresponding TDEP results, with significantly higher N_r deposition in the agricultural lands to the west of the GYA and hot spots due to wildfires that are not evident in the TDEP results. Both sets of results showed higher N_r wet deposition at the higher-elevation sites in the interior of the GYA, which was associated with higher precipitation rates. However, the TDEP N_r wet deposition was generally higher throughout the GYA, with an annual average N_r wet deposition rate of 2.0 N ha⁻¹ yr⁻¹ vs. 1.3 N ha⁻¹ yr⁻¹ from CAMx. Precipitation maps generated by WRF and PRISM across the GYA had similar spatial patterns, with hot spots located in high-elevation mountain ranges, though the WRF annual precipitation rates were on average 73 % higher than the PRISM estimates.

The annual N_r deposition budget and the annual precipitation rate from TDEP and the CAMx simulations for eight Class I areas over the GYA are compared in Fig. 5. The reported CAMx dry and wet N_r deposition values in Fig. 5 are the averages of the simulation values at corresponding grid cells for each area. Generally, results from the CAMx model agreed well with TDEP results in terms of replicating spatial gradients and ratios of oxidized vs. reduced N deposition. The TDEP 2011 annual N_r deposition at the GYA receptor sites was in the range of 2.8–5.4 kg N ha⁻¹ yr⁻¹, while the corresponding values for CAMx were 2.2–4.3 kg N ha⁻¹ yr⁻¹. Both results showed the west-to-east gradient (Prenni et al., 2014) with higher N_r deposition at the western side of the GYA and relatively low values in the Fitzpatrick Wilderness. Also, both models showed the importance of reduced N_r in the GYA, with a nearly 50 % or higher contribution to the total N_r deposition budget. However, the two models differed on the ratio of dry vs. wet N_r deposition, with CAMx simulating a higher fraction from dry N_r deposition than TDEP.

4 Source apportionment of N_r deposition over the GYA in 2011

The seasonal modeled N_r deposition budgets averaged over the GYA are presented in Fig. 6. As shown, the total N_r deposition rates peaked in the summer (1.12 kg N ha⁻¹ season⁻¹) with somewhat lower rates in the spring (0.91 kg N ha⁻¹ season⁻¹) and fall (0.81 kg N ha⁻¹ season⁻¹) and with winter rates (0.29 kg N ha⁻¹ season⁻¹) being about a factor of 3 smaller than in the other seasons. These patterns are similar to the measured and modeled data presented in Fig. 3. In total, the annual model N_r deposition was 3.13 kg N ha⁻¹ yr⁻¹, with wet deposition accounting for only ∼40 %. Reduced N compounds were the largest contributor, except in winter, which is consistent with past studies (Li et al., 2017). Contributions from organic N compounds are not measured in routine monitoring programs. Together they accounted for <10 % of the N_r deposition, suggesting a small but significant contribution. This is also less than has been measured in field studies conducted at GTNP (Benedict et al., 2013a; Prenni et al., 2014) and in RMNP (Benedict et al., 2013b), where the GrandTReNDS study showed on average a 8 %–18 % contribution from organic N to total N_r deposition budgets during the whole campaign period and up to 39 % in June at the NOAA Climate Station site (Fig. 7 in Benedict et al., 2013a).

Figure 6 Seasonal CAMx simulated N_r deposition budgets averaged over the GYA in 2011. The left axis is the relative contribution of different N_r species to seasonal N_r deposition while the right axis corresponds to the black diamonds for seasonal total N_r deposition as an absolute (kg N ha⁻¹).

Figure 7 Contributions of source sectors to the mean total N_r deposition, dry N_r deposition, and wet N_r deposition over the GYA in different seasons in 2011. (a) The source sector contributions as an absolute and (b) the corresponding contributions as a percentage.

Figure 8 Seasonal patterns of different source sectors' (agriculture, oil and gas activities, fires, others (e.g., anthropogenic, biogenic, lightning, and boundary conditions)) contributions to total N_r deposition over the GYA in 2011. The first column is the total seasonal N_r deposition patterns in kg N ha⁻¹ while the following five columns are the seasonal patterns of relative contributions from different source sectors.

The relative contributions from the four modeled source sectors (AG, OG, Fire, and Other) and the BCs averaged over the GYA are presented in Fig. 7, while Fig. 8 presents the seasonal and spatial patterns of their contributions over the GYA. As shown in Table S2, the AG source sector was composed of almost all reduced N compounds (>99 %), while the Other source sector was primarily composed (97 %) of oxidized N compounds, with about 88 % originating from anthropogenic combustion emissions, including point and mobile sources, and the remainder originating from the natural emissions from soil and lightning. Contributions from the Fire and the BCs sectors were more evenly split between reduced and oxidized N contributions.

Reduced N from the AG source sector was the largest contributor in the spring (40 %) and fall (41 %) seasons, while oxidized N from the Other source sector was the largest contributor in summer (29 %) and winter (44 %) (Fig. 7). In terms of geographic impact (Fig. 8), AG emissions contributed as much as 80 % of the total N_r deposition in the western portion of the GYA during the spring and fall, which was associated with the outflow from the Snake River valley. In the model, NH₃ from regional agriculture activities was treated as being from surface area sources (i.e., emitted into the first model layer, which is approximately 24 m thick). These low-level emissions can be quickly deposited to the surface unless there is sufficient vertical mixing to inject the NH₃ into the upper levels of the atmosphere (Ferm, 1998; Fenn et al., 2003) or if it reacted with acidic gases and aerosols. Consequently, it is likely that a higher fraction of the modeled NH₃ emissions from the AG sector will be deposited in the lower-elevation periphery of the GYA near the agricultural lands and not impact the more distant mountainous interior (Fig. 2). The incorporation of the bidirectional NH₃ flux could extend the NH₃ emission footprint (Bash et al., 2013; Zhu et al., 2015).

The OG source sector contributed only about 1 % of the total N_r deposition over the GYA, with contributions of 10 % or more occurring during winter in the southeastern corner of the GYA where nearby OG activity in the Jonah Field and Pinedale Anticline was taking place. Wildfires are episodic and their locations and magnitudes vary significantly from year to year (Westerling and Swetnam, 2003; Parisien et al., 2012). In 2011, fire events contributed on average 18 % of the total N_r deposition in the GYA. Most of the wildfire happened in summer and fall, while agriculture and prescribed burning occurred in winter and spring. Near the fire activities, the contribution to N_r deposition could be more than 90 %, as seen in Fig. 8. The footprint of fire emission impacts depends on the simulated injection height of the fire plumes. The emissions from fires that occurred within the GYA during the summer and fall likely remained within the mixed layer and had less of a chance to be transported far downwind to impact more distant areas (Fig. S4). The Other source sector had relatively uniform contributions throughout the GYA, indicative of contributions from regional sources. The Other sector accounted for 26 % of the annual N_r deposition, with its largest absolute contributions in the summer, but had the highest relative contribution in the winter at 44 % when AG contributions were at their lowest. Finally, the BCs had high contributions, often over 20 %, with the highest contributions occurring in the northern part of the GYA and at higher-elevation sites.

The seasonal contributions from the modeled source regions and sectors to the average total N_r deposition over the GYA are summarized in Fig. 9. As shown, the Snake River valley in Idaho was the largest contributor (in all seasons), with annual mean contributions of 38 % and a maximum contribution of 43 % in fall. Most (74 %) of the N_r from this region was from the AG source sector and was composed of reduced N (Table S4). The next four largest contributors, on average, were the BCs (21 %), western Wyoming (8 %), California (7 %), and northern Utah (6 %). The impact of emissions from Wyoming on the GYA during summer and fall (14 % and 16 %, respectively) was more pronounced than winter and spring (5 % and 7 %, respectively). The contributions of long-range transport from California and the BCs were higher during spring and winter.

Seasonal source apportionment results of the average dry and wet N_r deposition over the GYA are shown in Figs. 7 and 9. Compared to the results for total N_r deposition, the dry N_r deposition had higher contributions from closer sources, such as the Snake River valley (46 % for dry vs. 38 % for total), with emissions primarily from AG sources. Similarly, contributions to dry N_r deposition from Wyoming were 15 % compared to 12 % for total N_r deposition and ranked as the second-largest contributor. The contributions from distant source regions decreased. For example, the BCs decreased from 21 % for total N_r deposition to 12 % for dry N_r deposition.

Figure 9 Contributions of source regions to the mean total N_r deposition, dry N_r deposition, and wet N_r deposition over the GYA in different seasons in 2011. (a) The source region contributions as an absolute and (b) the corresponding contributions as a percentage.

Figure 10 Contributions of different source sectors as well as boundary conditions for total N_r deposition in 2011 at 10 points of interest for critical load exceedance (see Table 2 for site locations and ecosystem impacts). The black-and-white pies are the contributions by source sector while the color pies are the contributions by source region. The color contour for the GYA boundary is the terrain heights with the legend at at the right.

Table 2 Total reactive nitrogen (N_r) deposition and critical loads for receptor points in the Greater Yellowstone Area in Wyoming.

¹ Comprehensive Air Quality Model with extensions. ² NADP total deposition maps. ³ The range of critical loads (CLs) for different effects on the selected sensitive ecosystem receptor is from US CLAD (Critical Loads for Sulfur and Nitrogen Access Database), version 2.5 (Lynch et al., 2015). The level of confidence is based on the work of Pardo et al. (2011). The lower ends of the range were used in this study as a measured CL. ⁴ The CL values were for lichen response at sites with snowpack as a sensitive ecosystem.

The opposite pattern is seen for wet N_r deposition, where the contributions from the distant source regions increased relative to the neighboring ones. The annual contributions from the BCs increased to 34 % and peaked in spring and summer at 37 %, associated with higher precipitation amounts than the other two seasons. Annual contributions from sources in California (10 %) and Utah (8 %) surpassed Wyoming (7 %). Furthermore, the seasonal variation for wet N_r deposition was different from dry and total N_r deposition, with the highest deposition rates occurring in spring as opposed to summer.

The GYA has been the focus of several ecological assessments of the response of ecosystems to changing N_r deposition levels (Spaulding et al., 2015; Nanus et al., 2017). Figure 10 presents the source attribution results for 10 sites within the GYA where either ecosystem response studies or deposition monitoring has been conducted for lichen diversity, alpine lake chemistry, and snowpack analysis. In Table 2, the CL values are provided as a range of lower-end and upper-end estimates of the annual total inorganic N_r deposition values (Lynch et al., 2015) with confidence levels (Pardo et al., 2011). The simulated N_r deposition exceeded the lower CL values at three of the 10 sites, specifically, Holly Lake, Pinedale, and Tower Falls. Comparatively, the 2011 TDEP N_r deposition results exceeded the CL at 6 out of 10 sites (Black Joe Lake, Biscuit Basin, Holly Lake, Jedediah Smith Wilderness, Pinedale, and Twin Island). As shown in Fig. 10, the sites that exceeded the CL tend to be in high-alpine locations, with four of these sites on the western slope of the mountains, which are downwind of the Snake River valley. These results are consistent with another modeling study to access CL exceedances in Class I areas using GEOS-Chem (Ellis et al., 2013; Lee et al., 2016). In addition, in one study (Nanus et al., 2017) over 30 % of the GYA was estimated to potentially exceed lower N_r deposition CL thresholds, with the greatest impacts in sensitive high-elevation basins, including areas within national parks and wildernesses.

In terms of emission sectors and source regions contributing to the total annual N_r deposition at CL exceedance sites, emission sources from the Snake River valley were the largest contributors (27 %–32 %), and AG emissions were the largest source of this subset. The next three largest contributors were transport from the BCs (23 %–25 %) and emissions from northern Utah (8 %–15 %) and California (7 %–8 %). Wyoming emissions associated with the OG and Fire emission sectors contributed around 3 %–5 % and 14 %–23 %, respectively, of the N_r budget for receptor sites at the southeastern corner of the GYA.

5 The influence of model bias on source apportionment results

It is evident from the results in Sect. 4 that the attribution of total N_r deposition to source regions and sectors is sensitive to NH₃ dry deposition rates, the relative contributions of dry and wet deposition, and the concentrations of N compounds from the BCs. However, the model evaluation revealed a significant underestimation of NH₃ concentrations and overestimation of HNO₃ concentrations and precipitation rates; thus, these modeling errors could bias the source attribution results. To better understand the potential effects of these biases, sensitivity analyses of the source attributions to changes in NH₃ dry deposition rates and average precipitation rates as well as potential biases in the BCs were evaluated.

To test the sensitivity of the apportionment to NH₃ dry deposition rates, the deposition velocities were reduced by increasing the NH₃ resistance scaling factor by 10 %, following the methodology used in Thompson et al. (2015). The Zhang et al. (2003) dry deposition scheme was used in the CAMx simulations (Table S1), and this resistance scaling factor is designed to address the rapid removal of “sticky” compounds such as HNO₃ and NH₃ and can yield a nonlinear response in the estimated dry deposition velocity. July and August 2011 were simulated using the modified deposition velocity, and these results will be referred to as “DV_0.1”. The 10 % change in the resistance factor slowed the NH₃ deposition velocity from 2.5–4 to 1–1.5 cm s⁻¹ over the GYA, resulting in values more comparable to those used in the GrandTReNDS study (Benedict et al., 2013a; Prenni et al., 2014). The simulated NH₃ concentrations for the DV_0.1 case increased throughout the GYA compared to the base case. This resulted in better agreement with NH₃ measurements at the Grand Targhee and NOAA climate station sites but poorer agreement at the Driggs monitoring site (Fig. S5). The slower dry deposition velocities result in a longer NH₃ lifetime, allowing it to travel farther from nearby source regions, e.g., the Snake River valley, into the GYA and cause a more homogeneous concentration pattern throughout the GYA (Fig. S6). As shown in Fig. 11, the slower deposition velocities also somewhat altered the source attribution results. The contribution from the AG emission sector increased with the DV_0.1 simulation to 23 % compared to 19 % in the base case, with a smaller decrease in the contributions from the Other and the Fire sectors. This change was due to small increases in the contributions from the Snake River valley and northern Utah and decreases from Wyoming. Overall, decreasing the NH₃ dry deposition rate by about a factor of 2 had only a small impact on the N_r deposition budget and source apportionment results in the GYA. It is important to note that, although this was a significant reduction in the simulated dry deposition velocity for NH₃, it still represents a relatively rapid removal rate compared to other species, and NH₃ is quickly lost from the atmosphere in either case. It is known that NH₃ deposition in many environments is a bidirectional as opposed to a unidirectional process, and modeling the NH₃ flux as a bidirectional process may further decrease the bias for ambient NH_x concentration simulations (Bash et al., 2013; Wen et al., 2014; Whaley et al., 2018). The key process in air quality models to represent the re-emission of NH₃ from soil and plants to the atmosphere is the estimation of the available soil NH_x pool and the parameterization of compensation points for the conditions to re-emit NH₃ (Zhang et al., 2010; Whaley et al., 2018). In the CMAQ model, the bidirectional NH₃ deposition was realized by coupling with the United States Department of Agriculture's (USDA) Environmental Policy Integrated Climate (EPIC) agroecosystem model to provide the fertilization timing, rate, and composition (Bash et al., 2013). There is no similar parameterization available in the current CAMx model. Furthermore, the CAMx source apportionment tools cannot properly account for the origin of NH₃ concentrations at a receptor that has been deposited and then re-emitted.

Figure 11 The sensitivity of NH₃ dry deposition velocity (a, b “base” case; c, d “DV_0.1” case with NH₃ dry deposition velocity slowing down) to source apportionment results over the GYA during July–August 2011. Panels (a) and (c) show the contributions by source region as an absolute and as a percentage while (b) and (d) show the contributions by source sector.

Figure 12 Ratio of simulated vs. measured particulate nitrate (PNO₃) concentrations against the boundary contributions to simulated PNO₃ at IMPROVE sites over a 12 km domain.

The CAMx simulation overestimated the wet N_r deposition at measured sites, which was likely associated with an overestimation in the precipitation rates from WRF, especially at high-elevation sites. This precipitation rate bias was large, with the annual precipitation over the GYA more than 73 % higher than the PRISM estimates. We used the Noah land-surface model and Kain–Fritsch scheme cumulus parametrization in the WRF simulations (Table S1), and those physical module configurations were reported to have the tendency to overestimate precipitation (Warrach-Sagi et al., 2013). To evaluate the impact of the overestimation in precipitation on the source attribution results, the seasonal wet deposition rates were scaled to the measured precipitation rates at all NADP NTN and GrandTReNDS monitoring sites, following the procedures by Appel et al. (2011). This was equivalent to scaling the modeled wet deposition rates by the ratio of the measured to modeled precipitation rates. This approach assumes that the concentrations of N_r in the precipitation were the same in the model and measured data, which was not the case. After the precipitation adjustment, the correlation between the simulated and measured N_r wet deposition improved (Fig. S7). Within the GYA, however, the scaled N_r wet deposition underestimated the measured N deposition by about a factor of 2 and significantly underestimated the ratio of wet to dry deposition. Consequently, scaled wet deposition results were not used in this assessment. The overestimation of HNO₃ concentrations in the GYA is another reason for the wet N_r deposition overestimation. However, its impact on source apportionment results was not conducted here due to unclear reasons for the model bias (emission, chemistry, meteorology, deposition scheme) and limited computational resources.

The BC used in this work was derived from a MOZART global model simulation. An alternative set of BCs from the GEOS-Chem global model was also evaluated. Both sets of BCs resulted in high contributions to the total N_r deposition in the GYA, with the GEOS-Chem results having a slightly higher average contribution of 23 % compared to 21 % for MOZART (Figure S8). However, the GEOS-Chem BCs resulted in higher relative contributions of oxidized N to the total N_r deposition rate compared to the MOZART BCs (51 % and 45 %, respectively). The poor correspondence in the oxidized to reduced N_r split is reflective of the large uncertainties in the BC contributions to the N_r deposition and suggests that more evaluation of the global model results is warranted.

To examine the potential bias in the BC contributions, the simulated PNO₃ concentrations were compared to measurements from the IMPROVE monitoring program over the western United States for 2011. This comparison is shown in Fig. 12, in which the ratio of the simulated to measured PNO₃, i.e., an estimate of the bias, is plotted against the relative fraction of the contribution of the BCs to the simulated PNO₃. The data were first segregated by the fractional contribution of the BCs and then averaged together. As shown, for the MOZART BCs, the bias increased with larger relative contributions from the BCs, and when the BC fraction was 60 %, the bias was more than a factor of 2. This suggests that at least the particulate nitrate concentrations from the BCs are overestimated and possibly other N_r compounds from the BCs as well. In a CMAQ simulation using BCs derived from a GEOS-Chem simulation, Baker et al. (2015) also found that the contributions from the BCs to PNO₃ were overestimated when compared to IMPROVE data.

6 Summary and discussion

The CAMx model and its PSAT source apportionment tool were used to examine and quantify the contributions of four different source sectors and 27 source regions and the boundary conditions (BCs) to the 2011 total inorganic N_r deposition within the GYA. The source sectors were agriculture (AG), oil and gas activities (OG), wild and prescribed fires (Fire), and remaining contributions labeled as Other. The Other sector was primarily composed of oxidized N originating from anthropogenic combustion sources, including mobile and point sources, and the AG sector was almost entirely composed of reduced N compounds. Fire and the BCs were a mix of reduced and oxidized N compounds. This assessment focused on only the inorganic N fraction. There is measured evidence that organic N (Benedict et al., 2013a; Prenni et al., 2014) is a significant contributor to N_r deposition, and the inability to assess its origin in the current CTM is an important uncertainty in this work. Nevertheless, this N_r source apportionment work is the first thorough analysis of the origin of inorganic N_r in the GYA using a regional air quality modeling platform. The detailed source sector and source region configurations in PSAT enabled quantitative, though uncertain, estimates of their relative importance. This information is needed by stakeholders and regulators to understand the causes of excess N_r deposition in the GYA, monitor changes in N_r deposition, and develop possible future mitigation strategies.

Overall, the model simulation had a reasonable capacity to reproduce the measured seasonal and annual total N_r deposition levels throughout the GYA. However, the model simulation underestimated the measured NH₃ concentrations by 65 % on average and overestimated the measured HNO₃ by 108 %. Therefore the model tended to overestimate the contributions from oxidized N compounds and underestimate those from reduced N compounds to total N_r deposition. In addition, both reduced and oxidized N_r wet depositions were overestimated by 20 %–30 %, which was due, at least partially, to the simulated precipitation frequency and magnitude being too high in the model. These biases suggest that the modeled contributions from the AG emission sector were underestimated, while those from the Other sector's activities were overestimated.

The simulated annual total N_r deposition over the GYA in 2011 was 3.13 kg N ha⁻¹ yr⁻¹ and exceeded the CL estimates for lichen and lake chemistry primarily at high-elevation sites on the western slope and southern portion of the GYA. This finding is consistent with other studies using global models. Ellis et al. (2013) used the GEOS-Chem model to estimate the N_r deposition to Class I areas for 2006 and showed that the simulated total N_r deposition at GTNP (2.9 kg N ha⁻¹ yr⁻¹) and YNP (2.6 kg N ha⁻¹ yr⁻¹) exceeded the low end of CL for lichens (2.5 kg N ha⁻¹ yr⁻¹).

Emissions from the AG sector within the modeling domain were the largest contributor to the GYA total N_r deposition budget at 34 % year⁻¹. The contributions from the Other sector were also large at 26 %. The OG emission sector generally had a small contribution, except at the southern edge of the GYA, where it could contribute over 10 % of the total N_r deposition during winter months, with almost half of the OG contributions originating from emissions in the neighboring Jonah Field in western Wyoming. The Fire emission sector also had a significant contribution of 18 % over the year. This was due to regional contributions from fires throughout the west and large contributions (>90 %) in areas within the GYA where several wildfires occurred (Fig. 8). The large impact from fires within the GYA is notable since the episodic nature of fire will result in differing year-to-year contributions from this uncontrollable sector.

The largest impact from the AG emission sector originated from sources relatively close to the GYA, and the Snake River valley accounted for 74 % of the annual agricultural contribution. The agricultural contribution from Wyoming was 7 %, and more-distant source regions in northern Utah, California, and the northwestern United States each accounted for 4 %–5 % of the agricultural contribution. Nearly half (45 %) of the N_r deposition from the OG emission sector originated within Wyoming, especially the Upper Green River (27 %). The largest impact from the Fire emission sector originated from the Snake River valley (33 %) and within the GYA (25 %). The Other emission sector was more evenly distributed among near and distant regions, with the Snake River valley 23 %, Wyoming 17 %, and northern Utah, California, and the northwestern United States accounting for 14 %–16 % of the N_r deposition.

Long-range transport of N species from the BCs, which primarily originated from international sources, contributed 21 % of the total N_r deposition within the GYA during 2011 and had the largest absolute contribution during the summer. Several studies have shown the importance of international source contributions to particulates and N deposition within the continental United States (Park et al., 2004; Brewer and Moore, 2009; Zhang et al., 2012; Fann et al., 2013; Baker et al., 2015; Thompson et al., 2015). However, the BCs contribution in this work is on the high end of the reported values. For example, in a similar modeling study by Thompson et al. (2015), the estimated contribution of BCs to N_r deposition in Rocky Mountain National Park in 2009 was 13 %. Zhang et al. (2012) used the GEOS-Chem model to evaluate N deposition in the United States during 2006–2008 and showed that foreign anthropogenic contributions were generally <10 % but could rise up to 30 % near the Canadian and Mexican borders. In addition, our evaluations of the BCs suggest that the contribution of the BCs to ambient PNO₃ and possibly other N_r compounds was overestimated (Fig. 12), clearly suggesting that more research is needed on the role of distant emission sources in impacting N deposition in remote areas, as well as further investigations into model biases.

The observed precipitation in 2011 was ∼30 %–50 % higher than the historical average (NOAA, 2012), with the largest bias occurring at the eastern sites in the GYA (Fig. S9). This suggests that dry deposition of NH₃ may be a more important contributor to total N_r deposition during spring than that observed during GrandTReNDS. Also, considering that the wet deposition in the GYA tended to be overestimated and the precipitation amount in 2011 was anomalously high, the source regions identified as having a higher weighting on the annual wet N_r deposition budget (e.g., California) may not have such a significant impact as the current PSAT results suggested.

As discussed, source apportionment assessments of N_r and its deposition to remote, ecologically sensitive areas such as the GYA have large uncertainties. Many of these uncertainties are known to the air quality modeling community, including the challenges of simulating precipitation in complex terrain, adequately characterizing NH₃ emissions from agricultural operations, the occurrence of wildfires, and the difficulty in simulating the NH₃ bidirectional flux and the deposition flux of the other N_r compounds. Contributions from long-range transport of international emissions can also play a significant role in deposition in remote locations in the western United States. Further refinement in all of these areas is required to better understand and estimate the relative contributions of emission sources to excess N deposition within the GYA. Nevertheless, the modeling assessment showed that reduced N contributed more than 50 % of the total N_r deposition over the GYA, with >90 % of the NH₃ emissions originating from agriculture sources. In addition, the Snake River valley in Idaho accounted for 74 % of the agricultural contribution to the total N_r deposition. Significant contributions from more distant sources, e.g., California and international sources, to both oxidized and reduced N_r deposition illustrate the regional nature of the N_r deposition problem. Emissions of oxidized N compounds are projected to continue to decrease, while emissions of ammonia are projected to remain relatively constant or increase (Li et al., 2016). This will further increase the importance of the AG sector. However, exceedances of CL are still relatively small, and it is possible that decreased oxidized N deposition could reduce the N_r deposition sufficiently to bring total N_r deposition below the CL in some GYA ecosystems.