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

Abstract

Research has shown that excess reactive nitrogen (N${}_{\text{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${}_{\text{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${}_{\text{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${}_{\text{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${}_{\mathrm{3}})$ 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${}_{\text{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${}_{\text{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.

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${}_{\text{r}}$) compounds in air, rain, and snowpack samples over the GYA have been reported in the past 30 years and linked to N${}_{\text{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${}_{\text{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${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$, compared with 1.5–2.5 kg N ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$ in 2000 (NADP, 2016). This is relevant to the long-term conservation of the area because as N${}_{\text{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${}_{\text{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${}_{\text{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${}_{\text{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 (${\mathrm{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 ${\mathrm{NO}}_x$ result in nitric acid (${\mathrm{HNO}}_{\mathrm{3}}$), particulate nitrate (${\mathrm{PNO}}_{\mathrm{3}}$), and other compounds. Reduced N arises primarily from ammonia (${\mathrm{NH}}_{\mathrm{3}}$) 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 ${\mathrm{NH}}_{\mathrm{3}}$ 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 ${\mathrm{NH}}_{\mathrm{3}}$ 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 ${\mathrm{NO}}_x$ with volatile organic compounds (Cape et al., 2011; Reay et al., 2012). With the steady decline of ${\mathrm{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${}_{\text{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 ${\mathrm{NH}}_{\mathrm{3}}$ concentrations, with higher concentrations west of the Teton mountain range. Concurrently measured concentrations of ${\mathrm{HNO}}_{\mathrm{3}}$, PM${}_{\mathrm{2.5}}$ (particulate matter with an aerodynamic diameter of less than 2.5 $\mathrm{\mu }$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${}_{\text{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${}_{\text{r}}$ deposition over the Teton area. Lee et al. (2016) used the adjoint version of GEOS-Chem to quantify the sources of N${}_{\text{r}}$ deposition in eight selected federal Class I areas in 2010 and found a nonnegligible footprint ($>\mathrm{20}$ %) of N${}_{\text{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 ${\mathrm{NO}}_x$ and livestock ${\mathrm{NH}}_{\mathrm{3}}$ were also found to be major sources of N${}_{\text{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${}_{\text{r}}$ were larger than those from sources of oxidized N${}_{\text{r}}$ (Thompson et al., 2015; Malm et al., 2016).

In this work, we add to the growing body of N${}_{\text{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${}_{\text{r}}$ throughout the GYA. Compared with previous N${}_{\text{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${}_{\text{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${}_{\text{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${}_{\text{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.

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 omitted. See PDF]

Annual ${\mathrm{NH}}_{\mathrm{3}}$ 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 (${}^{\text{a}}$ Ammonia Monitoring Network; ${}^{\text{b}}$ Clean Air Status and Trends Network; ${}^{\text{c}}$ Grand Teton Reactive Nitrogen Deposition Study; ${}^{\text{d}}$ Interagency Monitoring of Protected Visual Environments; ${}^{\text{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.

[Figure omitted. See PDF]

Modeling system for N${}_{\text{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 ${\mathrm{NH}}_{\mathrm{3}}$ and ${\mathrm{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 ${\mathrm{NH}}_{\mathrm{3}}$ in the inner 12 km domain as well as the monitoring sites or receptor areas used for the model evaluation and analysis. For ${\mathrm{NH}}_{\mathrm{3}}$ 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 ${\mathrm{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${}_{\text{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${}_{\text{r}}$ deposition, but they can mask important contributions from other sources that are significant in nonfire years.

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.

						Site	N${}^{\mathrm{10}}$
Species		Network	Duration	OBS${}^{\mathrm{7}}$	SIM${}^{\mathrm{8}}$	no.${}^{\mathrm{9}}$	(% completeness)	$R^{\mathrm{11}}$	NMB${}^{\mathrm{12}}$	NME${}^{\mathrm{13}}$	FB${}^{\mathrm{14}}$	FE${}^{\mathrm{15}}$
Concentration	${\mathrm{NH}}_{\mathrm{3}}$ (ppb)	AMoN${}^{\mathrm{1}}$	22 Sep–12 Dec	0.49	0.30	1	7 (100 %)	0.20	$-\mathrm{65}$ %	67 %	$-\mathrm{52}$ %	53 %
		GrandTReNDS${}^{\mathrm{2}}$	5 Apr–21 Sep	0.55	0.46	3	434 (97.7 %)	0.30	$-\mathrm{16}$ %	57 %	$-\mathrm{42}$ %	63 %
	${\mathrm{HNO}}_{\mathrm{3}}$ (ppb)	CASTNet${}^{\mathrm{3}}$	4 Jan–27 Dec	0.23	0.47	2	83 (98.8 %)	0.72	108 %	117 %	60 %	71 %
		GrandTReNDS${}^{\mathrm{2}}$	5 Apr–21 Sep	0.28	0.54	3	435 (97.9 %)	0.60	106 %	109 %	63 %	68 %
	${\mathrm{PNO}}_{\mathrm{3}}$ ($\mathrm{\mu }$g m${}^{-\mathrm{3}})$	CASTNet${}^{\mathrm{3}}$	4 Jan–27 Dec	0.19	0.25	2	83 (98.8 %)	0.42	37 %	76 %	26 %	64 %
		IMPROVE${}^{\mathrm{4}}$	3 Jan–29 Dec	0.14	0.22	4	332 (68.5 %)	0.35	58 %	108 %	51 %	80 %
		GrandTReNDS${}^{\mathrm{2}}$	5 Apr–21 Sep	0.13	0.15	3	435 (97.9 %)	0.45	15 %	71 %	14 %	60 %
	${\mathrm{PNH}}_{\mathrm{4}}$ ($\mathrm{\mu }$g m${}^{-\mathrm{3}})$	CASTNet${}^{\mathrm{3}}$	4 Jan–27 Dec	0.17	0.18	2	83 (98.8 %)	0.28	3 %	39 %	7 %	41 %
		GrandTReNDS${}^{\mathrm{2}}$	5 Apr–21 Sep	0.14	0.17	3	433 (97.7 %)	0.12	23 %	64 %	34 %	61 %
	${\mathrm{NH}}_x$ ($\mathrm{\mu }$g m${}^{-\mathrm{3}}{)}^{\mathrm{5}}$	GrandTReNDS${}^{\mathrm{2}}$	5 Apr–21 Sep	0.68	0.63	3	427 (96.2 %)	0.26	$-\mathrm{7}$ %	48 %	$-\mathrm{22}$ %	46 %
N deposition	${\mathrm{HNO}}_{\mathrm{3}}$ dry (kg N ha${}^{-\mathrm{1}})$	CASTNet${}^{\mathrm{3}}$	4 Jan–27 Dec	0.071	0.187	2	83 (98.8 %)	0.81	153 %	156 %	77 %	82 %
		GrandTReNDS${}^{\mathrm{2}}$	5 Apr–21 Sep	0.016	0.049	3	435 (97.9 %)	0.66	204 %	209 %	101 %	104 %
	${\mathrm{PNO}}_{\mathrm{3}}$ dry (kg N ha${}^{-\mathrm{1}})$	CASTNet${}^{\mathrm{3}}$	4 Jan–27 Dec	0.012	0.023	2	83 (98.8 %)	0.14	96 %	148 %	48 %	97 %
		GrandTReNDS${}^{\mathrm{2}}$	5 Apr–21 Sep	0.010	0.011	3	435 (97.9 %)	0.61	8 %	58 %	1 %	65 %
	${\mathrm{PNH}}_{\mathrm{4}}$ dry (kg N ha${}^{-\mathrm{1}})$	CASTNet${}^{\mathrm{3}}$	4 Jan–27 Dec	0.018	0.019	2	83 (98.8 %)	0.1	7 %	57 %	22 %	61 %
		GrandTReNDS${}^{\mathrm{2}}$	5 Apr–21 Sep	0.006	0.004	3	433 (97.7 %)	0.1	$-\mathrm{33}$ %	46 %	$-\mathrm{28}$ %	53 %
	NO3${}^{-}$ wet (kg N ha${}^{-\mathrm{1}})$	NTN${}^{\mathrm{6}}$	4 Jan–27 Dec	0.079	0.097	5	214 (82.3 %)	0.34	31 %	126 %	12 %	100 %
		GrandTReNDS${}^{\mathrm{2}}$	5 Apr–21 Sep	0.051	0.083	3	427 (96.2 %)	0.15	60 %	94 %	42 %	71 %
	NH${}_{\mathrm{4}}^{+}$ wet (kg N ha${}^{-\mathrm{1}})$	NTN${}^{\mathrm{6}}$	4 Jan–27 Dec	0.088	0.126	5	214 (82.3 %)	0.32	49 %	142 %	19 %	106 %
		GrandTReNDS${}^{\mathrm{2}}$	5 Apr–21 Sep	0.103	0.147	3	427 (96.2 %)	0.48	42 %	72 %	30 %	64 %
Precipitation (cm)		NTN${}^{\mathrm{6}}$	4 Jan–27 Dec	0.77	2.34	5	214 (82.3 %)	0.54	215 %	242 %	64 %	118 %
		GrandTReNDS${}^{\mathrm{2}}$	5 Apr–21 Sep	0.33	0.95	3	427 (96.2 %)	0.42	187 %	207 %	69 %	94 %

${}^{\mathrm{1}}$ AMoN samples are collected for 2 weeks. ${}^{\mathrm{2}}$ CASTNet samples are collected for 1 week. ${}^{\mathrm{3}}$ GrandTReNDS samples are collected daily. ${}^{\mathrm{4}}$ IMPROVE 24 h samples are collected every 3 days. ${}^{\mathrm{5}}$ NH${}_x=$ ${\mathrm{NH}}_{\mathrm{3}}+{\mathrm{PNH}}_{\mathrm{4}}$. ${}^{\mathrm{6}}$ NTN samples are collected for 1 week. ${}^{\mathrm{7}}$ Average observation. ${}^{\mathrm{8}}$ Average simulation. ${}^{\mathrm{9}}$ Number of sites. ${}^{\mathrm{10}}$ Number of samples; the values in the parentheses are the percentage of valid samples used for model performance evaluation. ${}^{\mathrm{11}}$ Pearson's correlation coefficient. ${}^{\mathrm{12}}$ Normalized mean bias. ${}^{\mathrm{13}}$ Normalized mean error. ${}^{\mathrm{14}}$ Fractional bias. ${}^{\mathrm{15}}$ fractional errors.

CAMx–PSAT treats nitrogen-containing compounds as one of seven species: gaseous ${\mathrm{NH}}_{\mathrm{3}}$; particulate ammonium (${\mathrm{PNH}}_{\mathrm{4}}$); reactive gaseous nitrogen (RGN), which includes primary emissions of ${\mathrm{NO}}_x$, nitrous acid (HONO), nitrate radical (${\mathrm{NO}}_{\mathrm{3}}$), and dinitrogen pentoxide (${\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}$); gaseous nitric acid (${\mathrm{HNO}}_{\mathrm{3}}$); gaseous peroxy nitrogen (TPN), including peroxyacetyl nitrate (PAN) and peroxynitric acid (PNA); gas-phase organic nitrate (NTR); and particulate nitrate (${\mathrm{PNO}}_{\mathrm{3}}$). 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${}_{\text{r}}$ deposition as the sum of ${\mathrm{NH}}_{\mathrm{3}}$ and ${\mathrm{PNH}}_{\mathrm{4}}$ and the oxidized N${}_{\text{r}}$ deposition as the sum of RGN, ${\mathrm{HNO}}_{\mathrm{3}}$, ${\mathrm{PNO}}_{\mathrm{3}}$, TPN, and NTR in units of kg N ha${}^{-\mathrm{1}}$.

Evaluation of CAMx-simulated N${}_{\text{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${}_{\mathrm{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 ${\mathrm{NH}}_{\mathrm{3}}$ concentrations from the Ammonia Monitoring Network (AMoN) (http://nadp.sws.uiuc.edu/AMoN/, last access: 30 August 2018); nitric acid (${\mathrm{HNO}}_{\mathrm{3}}$), ${\mathrm{PNO}}_{\mathrm{3}}$, and ${\mathrm{PNH}}_{\mathrm{4}}$ 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); ${\mathrm{PNO}}_{\mathrm{3}}$ and ${\mathrm{PNH}}_{\mathrm{4}}$ concentrations from the Chemical Speciation Network (CSN) (https://www3.epa.gov/ttnamti1/speciepg.html, last access: 30 August 2018); ${\mathrm{PNO}}_{\mathrm{3}}$ 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).

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 ${\mathrm{HNO}}_{\mathrm{3}}$ with a normalized mean bias (NMB) of 108 % and significantly underestimated ${\mathrm{NH}}_{\mathrm{3}}$ concentrations with a NMB $=-\mathrm{65}$ %. While the model had skill in reproducing the daily variation in ${\mathrm{HNO}}_{\mathrm{3}}$, with a correlation coefficient of $r=\mathrm{0.71}$, it had little skill for ${\mathrm{NH}}_{\mathrm{3}}$, with $r=\mathrm{0.2}$. The overestimation of ${\mathrm{HNO}}_{\mathrm{3}}$ 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 ${\mathrm{HNO}}_{\mathrm{3}}$ may be due to the uncertainty for the ${\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}$ uptake coefficient setting for heterogeneous reactions (Foley et al., 2010). The poor ${\mathrm{NH}}_{\mathrm{3}}$ results may be related to the high uncertainty in the ${\mathrm{NH}}_{\mathrm{3}}$ emission inventory (Clarisse et al., 2009) and important missing physical mechanisms in the model, including the lack of bidirectional ${\mathrm{NH}}_{\mathrm{3}}$ 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 ${\mathrm{NH}}_{\mathrm{3}}$ flux mechanism in the model should increase ambient ${\mathrm{NH}}_{\mathrm{3}}$ concentrations in the GYA and thus decrease the large model underestimation of ${\mathrm{NH}}_{\mathrm{3}}$ concentrations.

For ${\mathrm{PNO}}_{\mathrm{3}}$ and ${\mathrm{PNH}}_{\mathrm{4}}$ simulations in the GYA, CAMx overestimated both species, with better performance for ${\mathrm{PNH}}_{\mathrm{4}}$ than ${\mathrm{PNO}}_{\mathrm{3}}$ (3 % vs. 37 %, respectively, in terms of NMB) and better agreement for ${\mathrm{PNO}}_{\mathrm{3}}$ 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\sim \mathrm{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.

Model performance for (a, b) seasonal average N${}_{\text{r}}$ concentration, (c, d) seasonal accumulated N${}_{\text{r}}$ deposition budget, and (e, f) seasonal accumulated precipitation amount at collocated location sites (YNP and Pinedale) over the GYA in 2011. ${}^{\mathrm{1}}$ Clean Air Status and Trends Network. ${}^{\mathrm{2}}$ Comprehensive Air Quality Model with extensions. ${}^{\mathrm{3}}$ National Trends Network. ${}^{\mathrm{4}}$ Parameter-elevation Regressions on Independent Slopes Model. ${}^{\mathrm{5}}$ Weather Research and Forecasting model.

[Figure omitted. See PDF]

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). ${}^{\mathrm{1}}$ Grand Teton Reactive Nitrogen Deposition Study. ${}^{\mathrm{2}}$ Comprehensive Air Quality Model with extensions. ${}^{\mathrm{3}}$ Weather Research and Forecasting model.

[Figure omitted. See PDF]

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 ${\mathrm{HNO}}_{\mathrm{3}}$ is evident in all seasons. Also evident is the poor simulation of the seasonality in N${}_{\text{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 ${\mathrm{HNO}}_{\mathrm{3}}$ and underestimation of ${\mathrm{NH}}_{\mathrm{3}}$. Resolution of these biases requires additional research and these biases need to be taken into account when interpreting the source attribution of N${}_{\text{r}}$ deposition within the GYA.

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${}_{\text{r}}$ deposition values at the eight Class I areas are the average of the simulation values at corresponding grid cells for each area.

[Figure omitted. See PDF]

Evaluation against GrandTReNDS data

The GrandTReNDS campaign provides a unique opportunity to evaluate the capability of CAMx to simulate the N${}_{\text{r}}$ compounds and deposition budget. Detailed measurements, including ${\mathrm{NH}}_{\mathrm{3}}$, were made at three sites that crossed GTNP from west to east: Driggs, in the foothills just west of GTNP (43.74${}^{\circ }$ N, $-\mathrm{111.87}$${}^{\circ }$ W, elevation 1947 m); Grand Targhee, an upper-elevation site on the western edge of GTNP (43.78${}^{\circ }$ N, $-\mathrm{110.94}$${}^{\circ }$ W, elevation 2722 m); and the National Oceanic and Atmospheric Administration (NOAA) climate station site on the eastern edge of GTNP (43.66${}^{\circ }$ N, $-\mathrm{110.71}$${}^{\circ }$ 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${}_{\text{r}}$ deposition rates both in the month-to-month variation and across the sites. The difference in the dry ${\mathrm{NH}}_{\mathrm{3}}$ 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 ${\mathrm{HNO}}_{\mathrm{3}}$ dry deposition rates, which were 2–3 times larger than those derived from the measured data. This is consistent with the systematic overestimation of ${\mathrm{HNO}}_{\mathrm{3}}$ concentrations (NMB $=$ 106 % in Table 1). Other biases also exist, including an underestimation in the ${\mathrm{NH}}_{\mathrm{3}}$ dry deposition, which was somewhat balanced by an overestimation in the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ wet deposition (NMB $=$ 60 %). The underestimation of ${\mathrm{NH}}_{\mathrm{3}}$ concentration still existed (NMB $=$ $-\mathrm{16}$ %), and one of the possible reasons may be due to the overestimation of ${\mathrm{HNO}}_{\mathrm{3}}$ in the model pushing excessive partitioning of ${\mathrm{NH}}_{\mathrm{3}}$ into the particle phase, which can be shown by the better model performance for ${\mathrm{NH}}_x$ (${\mathrm{NH}}_x$ $=$ ${\mathrm{NH}}_{\mathrm{3}}+{\mathrm{PNH}}_{\mathrm{4}}$) 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.

Evaluation against NADP TDEP

TDEP maps (Schwede and Lear, 2014) are widely used in the land management community to assess total N${}_{\text{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${}_{\text{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${}_{\text{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${}_{\text{r}}$ deposition in the Snake River valley, northern Utah, and across the Wyoming state border area near GTNP, with values $>\mathrm{5}$ kg N ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$. Within the GYA (Fig. S3), the CAMx simulation had higher dry N${}_{\text{r}}$ deposition, which was more spatially heterogeneous than the corresponding TDEP results, with significantly higher N${}_{\text{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${}_{\text{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${}_{\text{r}}$ wet deposition was generally higher throughout the GYA, with an annual average N${}_{\text{r}}$ wet deposition rate of 2.0 N ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$ vs. 1.3 N ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$ 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${}_{\text{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${}_{\text{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${}_{\text{r}}$ deposition at the GYA receptor sites was in the range of 2.8–5.4 kg N ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$, while the corresponding values for CAMx were 2.2–4.3 kg N ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$. Both results showed the west-to-east gradient (Prenni et al., 2014) with higher N${}_{\text{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${}_{\text{r}}$ in the GYA, with a nearly 50 % or higher contribution to the total N${}_{\text{r}}$ deposition budget. However, the two models differed on the ratio of dry vs. wet N${}_{\text{r}}$ deposition, with CAMx simulating a higher fraction from dry N${}_{\text{r}}$ deposition than TDEP.

Source apportionment of N${}_{\text{r}}$ deposition over the GYA in 2011

The seasonal modeled N${}_{\text{r}}$ deposition budgets averaged over the GYA are presented in Fig. 6. As shown, the total N${}_{\text{r}}$ deposition rates peaked in the summer (1.12 kg N ha${}^{-\mathrm{1}}$ season${}^{-\mathrm{1}})$ with somewhat lower rates in the spring (0.91 kg N ha${}^{-\mathrm{1}}$ season${}^{-\mathrm{1}})$ and fall (0.81 kg N ha${}^{-\mathrm{1}}$ season${}^{-\mathrm{1}})$ and with winter rates (0.29 kg N ha${}^{-\mathrm{1}}$ season${}^{-\mathrm{1}})$ 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${}_{\text{r}}$ deposition was 3.13 kg N ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$, with wet deposition accounting for only $\sim \mathrm{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 $<\mathrm{10}$ % of the N${}_{\text{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${}_{\text{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).

Seasonal CAMx simulated N${}_{\text{r}}$ deposition budgets averaged over the GYA in 2011. The left axis is the relative contribution of different N${}_{\text{r}}$ species to seasonal N${}_{\text{r}}$ deposition while the right axis corresponds to the black diamonds for seasonal total N${}_{\text{r}}$ deposition as an absolute (kg N ha${}^{-\mathrm{1}}$).

[Figure omitted. See PDF]

Contributions of source sectors to the mean total N${}_{\text{r}}$ deposition, dry N${}_{\text{r}}$ deposition, and wet N${}_{\text{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 omitted. See PDF]

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${}_{\text{r}}$ deposition over the GYA in 2011. The first column is the total seasonal N${}_{\text{r}}$ deposition patterns in kg N ha${}^{-\mathrm{1}}$ while the following five columns are the seasonal patterns of relative contributions from different source sectors.

[Figure omitted. See PDF]

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 ($>\mathrm{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${}_{\text{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, ${\mathrm{NH}}_{\mathrm{3}}$ 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 ${\mathrm{NH}}_{\mathrm{3}}$ 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 ${\mathrm{NH}}_{\mathrm{3}}$ 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 ${\mathrm{NH}}_{\mathrm{3}}$ flux could extend the ${\mathrm{NH}}_{\mathrm{3}}$ emission footprint (Bash et al., 2013; Zhu et al., 2015).

The OG source sector contributed only about 1 % of the total N${}_{\text{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${}_{\text{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${}_{\text{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${}_{\text{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${}_{\text{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${}_{\text{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${}_{\text{r}}$ deposition over the GYA are shown in Figs. 7 and 9. Compared to the results for total N${}_{\text{r}}$ deposition, the dry N${}_{\text{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${}_{\text{r}}$ deposition from Wyoming were 15 % compared to 12 % for total N${}_{\text{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${}_{\text{r}}$ deposition to 12 % for dry N${}_{\text{r}}$ deposition.

Contributions of source regions to the mean total N${}_{\text{r}}$ deposition, dry N${}_{\text{r}}$ deposition, and wet N${}_{\text{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 omitted. See PDF]

Contributions of different source sectors as well as boundary conditions for total N${}_{\text{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.

[Figure omitted. See PDF]

Total reactive nitrogen (N${}_{\text{r}})$ deposition and critical loads for receptor points in the Greater Yellowstone Area in Wyoming.

					Total N${}_{\text{r}}$
					deposition (kg N ha${}^{-\mathrm{1}})$		Critical load (kg N ha${}^{-\mathrm{1}}{)}^{\mathrm{3}}$
Site	Site name	Latitude/	Elevation	Sensitive	CAMx${}^{\mathrm{1}}$	TDEP${}^{\mathrm{2}}$	Range	Confidence
ID	(state)	longitude	(m)	ecosystem				level
1	Absaroka-Beartooth Wilderness (MT)	45.49${}^{\circ }$ N, 110.51${}^{\circ }$ W	2536	Lichen	1.93	2.80	3.02–4.89	Reliable
2	Twin Island(MT)	45.07${}^{\circ }$ N, 109.81${}^{\circ }$ W	2829	Lake chemistry	1.53	3.99	2.5–7.1	Fairly reliable
3	Tower Falls(WY)	44.92${}^{\circ }$ N, 110.42${}^{\circ }$ W	2457	Snowpack	3.8	1.87	2.93–4.81${}^{\mathrm{4}}$	Reliable
4	Moose Meadow(ID)	44.63${}^{\circ }$ N, 111.24${}^{\circ }$ W	1885	Snowpack	3.38	2.36	3.52–5.40${}^{\mathrm{4}}$	Reliable
5	Biscuit Basin(WY)	44.46${}^{\circ }$ N, 110.83${}^{\circ }$ W	2050	Snowpack	2.69	3.49	3.39–5.27${}^{\mathrm{4}}$	Reliable
6	Jedediah Smith Wilderness(WY)	43.79${}^{\circ }$ N, 110.94${}^{\circ }$ W	1944	Lichen	3.03	6.36	3.40–5.27	Reliable
7	Holly Lake(WY)	43.79${}^{\circ }$ N, 110.79${}^{\circ }$ W	2230	Lake chemistry	3.15	5.50	2.5–7.1	Fairly reliable
8	Fitzpatrick Wilderness(WY)	43.40${}^{\circ }$ N, 109.66${}^{\circ }$ W	2890	Lichen	1.79	1.86	3.41–5.29	Reliable
9	Pinedale(WY)	42.93${}^{\circ }$ N, 109.79${}^{\circ }$ W	2246	Lichen	3.39	2.67	2.66–4.53	Reliable
10	Black Joe Lake(WY)	42.74${}^{\circ }$ N, 109.16${}^{\circ }$ W	3133	Lake chemistry	2.32	3.56	2.5–7.1	Fairly reliable

${}^{\mathrm{1}}$ Comprehensive Air Quality Model with extensions. ${}^{\mathrm{2}}$ NADP total deposition maps. ${}^{\mathrm{3}}$ 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. ${}^{\mathrm{4}}$ The CL values were for lichen response at sites with snowpack as a sensitive ecosystem.

The opposite pattern is seen for wet N${}_{\text{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${}_{\text{r}}$ deposition was different from dry and total N${}_{\text{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${}_{\text{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${}_{\text{r}}$ deposition values (Lynch et al., 2015) with confidence levels (Pardo et al., 2011). The simulated N${}_{\text{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${}_{\text{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${}_{\text{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${}_{\text{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${}_{\text{r}}$ budget for receptor sites at the southeastern corner of the GYA.

The influence of model bias on source apportionment results

It is evident from the results in Sect. 4 that the attribution of total N${}_{\text{r}}$ deposition to source regions and sectors is sensitive to ${\mathrm{NH}}_{\mathrm{3}}$ 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 ${\mathrm{NH}}_{\mathrm{3}}$ concentrations and overestimation of ${\mathrm{HNO}}_{\mathrm{3}}$ 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 ${\mathrm{NH}}_{\mathrm{3}}$ 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 ${\mathrm{NH}}_{\mathrm{3}}$ dry deposition rates, the deposition velocities were reduced by increasing the ${\mathrm{NH}}_{\mathrm{3}}$ 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 ${\mathrm{HNO}}_{\mathrm{3}}$ and ${\mathrm{NH}}_{\mathrm{3}}$ 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 ${\mathrm{NH}}_{\mathrm{3}}$ deposition velocity from 2.5–4 to 1–1.5 cm s${}^{-\mathrm{1}}$ 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 ${\mathrm{NH}}_{\mathrm{3}}$ concentrations for the DV_0.1 case increased throughout the GYA compared to the base case. This resulted in better agreement with ${\mathrm{NH}}_{\mathrm{3}}$ 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 ${\mathrm{NH}}_{\mathrm{3}}$ 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 ${\mathrm{NH}}_{\mathrm{3}}$ dry deposition rate by about a factor of 2 had only a small impact on the N${}_{\text{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 ${\mathrm{NH}}_{\mathrm{3}}$, it still represents a relatively rapid removal rate compared to other species, and ${\mathrm{NH}}_{\mathrm{3}}$ is quickly lost from the atmosphere in either case. It is known that ${\mathrm{NH}}_{\mathrm{3}}$ deposition in many environments is a bidirectional as opposed to a unidirectional process, and modeling the ${\mathrm{NH}}_{\mathrm{3}}$ flux as a bidirectional process may further decrease the bias for ambient ${\mathrm{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 ${\mathrm{NH}}_{\mathrm{3}}$ from soil and plants to the atmosphere is the estimation of the available soil ${\mathrm{NH}}_x$ pool and the parameterization of compensation points for the conditions to re-emit ${\mathrm{NH}}_{\mathrm{3}}$ (Zhang et al., 2010; Whaley et al., 2018). In the CMAQ model, the bidirectional ${\mathrm{NH}}_{\mathrm{3}}$ 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 ${\mathrm{NH}}_{\mathrm{3}}$ concentrations at a receptor that has been deposited and then re-emitted.

The sensitivity of ${\mathrm{NH}}_{\mathrm{3}}$ dry deposition velocity (a, b “base” case; c, d “DV_0.1” case with ${\mathrm{NH}}_{\mathrm{3}}$ 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 omitted. See PDF]

Ratio of simulated vs. measured particulate nitrate (${\mathrm{PNO}}_{\mathrm{3}}$) concentrations against the boundary contributions to simulated ${\mathrm{PNO}}_{\mathrm{3}}$ at IMPROVE sites over a 12 km domain.

[Figure omitted. See PDF]

The CAMx simulation overestimated the wet N${}_{\text{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${}_{\text{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${}_{\text{r}}$ wet deposition improved (Fig. S7). Within the GYA, however, the scaled N${}_{\text{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 ${\mathrm{HNO}}_{\mathrm{3}}$ concentrations in the GYA is another reason for the wet N${}_{\text{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${}_{\text{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${}_{\text{r}}$ deposition rate compared to the MOZART BCs (51 % and 45 %, respectively). The poor correspondence in the oxidized to reduced N${}_{\text{r}}$ split is reflective of the large uncertainties in the BC contributions to the N${}_{\text{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 ${\mathrm{PNO}}_{\mathrm{3}}$ 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 ${\mathrm{PNO}}_{\mathrm{3}}$, i.e., an estimate of the bias, is plotted against the relative fraction of the contribution of the BCs to the simulated ${\mathrm{PNO}}_{\mathrm{3}}$. 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${}_{\text{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 ${\mathrm{PNO}}_{\mathrm{3}}$ were overestimated when compared to IMPROVE data.

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${}_{\text{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${}_{\text{r}}$ deposition, and the inability to assess its origin in the current CTM is an important uncertainty in this work. Nevertheless, this N${}_{\text{r}}$ source apportionment work is the first thorough analysis of the origin of inorganic N${}_{\text{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${}_{\text{r}}$ deposition in the GYA, monitor changes in N${}_{\text{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${}_{\text{r}}$ deposition levels throughout the GYA. However, the model simulation underestimated the measured ${\mathrm{NH}}_{\mathrm{3}}$ concentrations by 65 % on average and overestimated the measured ${\mathrm{HNO}}_{\mathrm{3}}$ by 108 %. Therefore the model tended to overestimate the contributions from oxidized N compounds and underestimate those from reduced N compounds to total N${}_{\text{r}}$ deposition. In addition, both reduced and oxidized N${}_{\text{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${}_{\text{r}}$ deposition over the GYA in 2011 was 3.13 kg N ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$ 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${}_{\text{r}}$ deposition to Class I areas for 2006 and showed that the simulated total N${}_{\text{r}}$ deposition at GTNP (2.9 kg N ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$) and YNP (2.6 kg N ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$) exceeded the low end of CL for lichens (2.5 kg N ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$).

Emissions from the AG sector within the modeling domain were the largest contributor to the GYA total N${}_{\text{r}}$ deposition budget at 34 % year${}^{-\mathrm{1}}$. 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${}_{\text{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 ($>\mathrm{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${}_{\text{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${}_{\text{r}}$ deposition.

Long-range transport of N species from the BCs, which primarily originated from international sources, contributed 21 % of the total N${}_{\text{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${}_{\text{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 $<\mathrm{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 ${\mathrm{PNO}}_{\mathrm{3}}$ and possibly other N${}_{\text{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 $\sim \mathrm{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 ${\mathrm{NH}}_{\mathrm{3}}$ may be a more important contributor to total N${}_{\text{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${}_{\text{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${}_{\text{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 ${\mathrm{NH}}_{\mathrm{3}}$ emissions from agricultural operations, the occurrence of wildfires, and the difficulty in simulating the ${\mathrm{NH}}_{\mathrm{3}}$ bidirectional flux and the deposition flux of the other N${}_{\text{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${}_{\text{r}}$ deposition over the GYA, with $>\mathrm{90}$ % of the ${\mathrm{NH}}_{\mathrm{3}}$ emissions originating from agriculture sources. In addition, the Snake River valley in Idaho accounted for 74 % of the agricultural contribution to the total N${}_{\text{r}}$ deposition. Significant contributions from more distant sources, e.g., California and international sources, to both oxidized and reduced N${}_{\text{r}}$ deposition illustrate the regional nature of the N${}_{\text{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${}_{\text{r}}$ deposition sufficiently to bring total N${}_{\text{r}}$ deposition below the CL in some GYA ecosystems.

Model data have been made available through the Mountain Scholar repository service. The nitrogen deposition source apportionment results over the Greater Yellowstone Area can be accessed at https://hdl.handle.net/10217/191136 (Zhang et al., 2018). Contact the corresponding author for any additional data requests.

The supplement related to this article is available online at: https://doi.org/10.5194/acp-18-12991-2018-supplement.

RZ, TMT conducted the source apportionment simulations. RZ, TMT, MGB, JLH, BAS performed results analysis and interpretation. RZ, MGB, JLH, BAS prepared the manuscript. MDB contributed the critical load excess assessment. JAM and WCM were involved in the discussion of model uncertainties and edited the drafts of the paper.

The authors declare that they have no conflict of interest.

Acknowledgements

This work was funded by the National Park Service Air Resources Division under cooperative agreement P17AC00773. The assumptions, findings, conclusions, judgments, and views presented herein are those of the authors and should not be interpreted as necessarily representing the National Park Service policies. The Interagency Monitoring of Protected Visual Environments (IMPROVE) is a collaborative association of state, tribal, and federal agencies and international partners. The U.S. Environmental Protection Agency is the primary funding source, with contracting and research support from the National Park Service. The Air Quality Group at the University of California, Davis, is the central analytical laboratory, with ion analysis provided by the Research Triangle Institute and carbon analysis provided by the Desert Research Institute. We acknowledge the Total Deposition (TDEP) Science Committee of the National Atmospheric Deposition Program (NADP) for their role in making the TDEP data and maps available. Edited by: Leiming Zhang Reviewed by: two anonymous referees