Advances in Meteorology

Volume 2014 (2014), Article ID 414015, 19 pages

http://dx.doi.org/10.1155/2014/414015

Research Article

Meteorological and Back Trajectory Modeling for the Rocky Mountain Atmospheric Nitrogen and Sulfur Study II

Kristi A. Gebhart,¹ William C. Malm,² Marco A. Rodriguez,³ Michael G. Barna,¹ Bret A. Schichtel,¹ Katherine B. Benedict,⁴ Jeffrey L. Collett Jr.,⁴ and Christian M. Carrico³

¹Air Resources Division, National Park Service, Fort Collins, CO 80523, USA

²Cooperative Institute for Research in the Atmosphere, Colorado State University, Fort Collins, CO 80523, USA

³AECOM, Inc., Fort Collins, CO 80525, USA

⁴Department of Atmospheric Science, Colorado State University, Fort Collins, CO 80523, USA

Received 17 January 2014; Accepted 28 March 2014; Published 8 May 2014

Academic Editor: Florinda Artuso

Copyright © 2014 Kristi A. Gebhart et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The Rocky Mountain Atmospheric Nitrogen and Sulfur (RoMANS II) study with field operations during November 2008 through November 2009 was designed to evaluate the composition and sources of reactive nitrogen in Rocky Mountain National Park, Colorado, USA. As part of RoMANS II, a mesoscale meteorological model was utilized to provide input for back trajectory and chemical transport models. Evaluation of the model's ability to capture important transport patterns in this region of complex terrain is discussed. Previous source-receptor studies of nitrogen in this region are also reviewed. Finally, results of several back trajectory analyses for RoMANS II are presented. The trajectory mass balance (TrMB) model, a receptor-based linear regression technique, was used to estimate mean source attributions of airborne ammonia concentrations during RoMANS II. Though ammonia concentrations are usually higher when there is transport from the east, the TrMB model estimates that, on average, areas to the west contribute a larger mean fraction of the ammonia. Possible reasons for this are discussed and include the greater frequency of westerly versus easterly winds, the possibility that ammonia is transported long distances as ammonium nitrate, and the difficulty of correctly modeling the transport winds in this area.

1. Introduction

Atmospheric deposition of reactive nitrogen in Rocky Mountain National Park, Colorado (RMNP), and surrounding areas of the Rocky Mountains has been the subject of research for 30 years (e.g., [1]). During this time, scientists have endeavored to understand the levels and chemical composition of deposited nitrogen [2, 3], the spatial and temporal trends [4, 5], the chemical and physical mechanisms by which nitrogen enters aquatic and ecological systems [6, 7], effects on biota [8, 9] and the dominant sources [10], emission rates [8, 11], and geographical and meteorological conditions (e.g., [11–16]) that cause deposition in sensitive alpine areas. This paper focuses on the meteorological and source attribution aspects.

Concentration and deposition measurements for the Rocky Mountain Nitrogen and Sulfur II (RoMANS II) study were conducted during November 2008 through November 2009. Others [2, 3, 17] have described these measurements in detail. The year-long RoMANS II was a follow-up to RoMANS I [11, 15, 16, 18], in which data were collected at the same location for two 6-week periods during the spring and summer of 2006. RMNP is in North Central Colorado, straddling the Continental Divide. The most detailed measurements for both studies were at the “core site,” 40.38 deg N, 105.546 deg W, and 2750 meters above mean sea level (MSL), on the eastern slope of the park as shown in Figure 1. To the east, the terrain drops rapidly to the relatively flat high plains, where elevations are 1500–1800 meters MSL, while some peaks to the west reach 4300 meters MSL. The most densely populated region of Colorado, including the city of Denver (1610 meters MSL), is the Front Range urban corridor, which runs north-south through the center of the state from Southern Wyoming to Southern Colorado along the transition between foothills and plains. The predominant agricultural activities are in Northeastern Colorado, with many along the North Platte River valley, including some in and near the Front Range urban corridor.

(a)

(b)

Figure 1: Views of the RoMANS core site and meteorological measurement sites in Estes Park, CO (yellow pins), generated with Google Earth. Both views are looking north. Panel (a) shows RMNP boundaries in heavy green, Colorado county boundaries in pale green, and major highways in yellow. The North Platte River valley is the dark green area mostly along Highway 34. Panel (b) shows the orientation of the local valleys and surrounding terrain. The core site and the profiler are approximately 12 km apart. The profiler and the RAWS sites are about 3.4 km apart. On the larger map, Cheyenne, WY, and Denver, CO, are approximately 162 km apart.

The relatively small-scale spatial inhomogeneity arising from the transition between high plains, foothills, and alpine peaks presents a challenge for meteorological models. The prevailing wind direction is westerly in most locations, but easterly upslope winds occur due to either diurnal mountain valley circulations with winds near the surface blowing up valley during the day and reversing at night or synoptic weather patterns including low-pressure (counterclockwise circulation) to the south or west or high pressure (clockwise circulation) to the north [19].

1.1. Studies of Nitrogen Source Attribution in Northern Colorado

Estimates of source attribution of nitrogen at RMNP and surrounding areas have been of interest for approximately 30 years. The earliest studies based hypotheses about the origins of nitrogen on spatial and temporal patterns in measured concentrations, deposition, and meteorology. As computing power increased, these valuable analyses continued, while the use of air mass trajectories and, more recently, mesoscale chemical transport models, were added to the tools available. Since the 1980s, the population, energy use, and mix of sources in Colorado and the Western United States have evolved, so currently significant sources may be different than in previous decades.

Lewis et al. [1] examined deposition of major cations, anions, particulate matter, and hydrogen ions at 42 sites in Colorado from May 1982 to May 1983. They inferred sources by examining spatial gradients and prevailing wind patterns. The greatest acidity was along the Continental Divide of the northern and southern thirds of the state, which they stated was because strong acids were more effectively neutralized by carbonates at lower elevations. Dominant sources of oxidized nitrogen were hypothesized to be power plants to the west and urban areas to the east of the mountains. Reduced nitrogen compounds were not measured.

Parrish [12] measured airborne NO and NO₂ at a high elevation site during the 1980s. They stated that though the Front Range is usually downwind, emissions from there were the predominant sources responsible for enhanced acidic deposition. They also hypothesized that the mountain-valley flow might be a mechanism to transport anthropogenic air pollutants from the boundary layer to the free troposphere.

Baron and Denning [13] determined that during the 1980s, strong acid anions, acidity, ammonium, and high salt concentrations in precipitation originated east of RMNP and were transported up valley by funneling winds or convective instability and that these winds influenced a lower elevation (2500 meters) site throughout the year and a high elevation (3200 meters) site most strongly during the summer.

Neff [33] describes the regionally complex meteorology and air quality modeling associated with the 1987-1988 Denver Brown Cloud study. He discusses the influence of northeasterly winds bringing ammonia from rural areas and the mixing of those emissions with and from the urbanized Front Range. He modeled a few days with a chemical transport model, but the model failed to resolve important mesoscale meteorological structures, especially boundary layer and moist processes at the 8 km horizontal resolution. He suggested that data assimilation was necessary. He also pointed out that emissions may fluctuate in response to meteorological conditions in ways not accounted for in modeled emissions, for example, less auto traffic but more wood burning during snow fall.

Sievering et al. [34] concluded that during 1991–1994 emissions from the Front Range had little impact at a high altitude (3540 meters) site in the alpine tundra, but that diurnal upslope flows brought these emissions to a lower site (3020 meters) in the subalpine forest, especially during the summer and during midday. They argued that the higher site was usually impacted by westerly winds in the free troposphere. Isentropic trajectories indicated that the alpine site was most often influenced by sources in Southern California and Northwestern Mexico on days with high nitrogen concentrations. They also found that Front Range emissions can be mixed with emissions from the west. High sulfate and nitrate concentrations in precipitation at the subalpine site were highly correlated with trajectories from the west, including the Four Corners region, Southern California, and parts of Northern Mexico.

Losleben et al. [14] examined 1984–1997 wet deposition data at two sites, a high alpine site and one lower and closer to the Front Range. Precipitation at the higher site was more acidic and had higher conductivity than the lower site. They hypothesized that chemical loading from western sources might be greater than previously thought because westerly winds predominate at the higher site. They concluded that the Front Range was an important source of pollutants in precipitation especially during winter precipitation events, but that the much higher frequency of westerly flow meant that sources to the west were cumulatively important. The most- and least-acidic episodes were associated with southwesterly and northwesterly flows, respectively. They also discussed meteorological differences between winter and summer precipitation events at the two sites, the influence of terrain, and the importance of accounting for duration and amount of precipitation when estimating sources of wet-deposited pollutants.

Heuer et al. [35] examined spatial and seasonal patterns in wet deposition and snow pack data in Colorado during 1992–1997 to infer sources of nitrogen and sulfur compounds. They found that high elevation sites were influenced by sources both east and west of the divide, with sources to the west having greater influence during the winter, while those to the east have more influence during summer. They concluded that, overall, sources to the east have greater influence on precipitation chemistry in the Colorado Rockies because summertime concentrations were significantly higher than winter.

Ingersoll et al. [36] examined back trajectories associated with snow events on the east side of RMNP during the 1998-1999 winters. They observed that concentrations of nitrate and sulfate in snow samples were the lowest in storms from the north and east and elevated for air masses from the west. The highest concentrations were in storms from the south, including metropolitan Denver.

Sievering et al. [37] stated that the probability of upslope flow occurring on any single day at a subalpine forest site was 40–60% during 1999, with the peak in August and minimum in May.

In a 2003 review of nitrogen deposition in Colorado and Southern Wyoming, Burns [4] concluded that evidence up to that point supported a significant eastern source for atmospheric nitrogen deposition and that this was consistent with population growth and energy use trends.

As part of the RoMANS I study, Rodriguez et al. [11] described chemical transport modeling for two multiweek periods during spring and summer 2006. The model was unable to adequately reproduce observed concentrations of nitrogen compounds, probably mostly due to uncertainty in emissions and deposition velocities. They noted that upslope easterly flow on the east side of RMNP is more common during precipitation than during dry hours for all seasons and all hours of the day. They also noted that the vertical depth, but not the frequency, of easterly flow was greater during summer than during spring.

Also as part of RoMANS I, Gebhart et al. [16] employed a regression technique using back trajectories to apportion concentrations of nitrogen and sulfur species measured on the east side of RMNP to selected source regions. The highest measured concentrations were usually associated with transport from the east. However, the model estimated that the season-long mean attributions for all nitrogen species for both seasons were greater for regions to the west. Again, this was thought to be due to the much greater frequency of westerly winds. The receptor model performed much better for spring than summer, indicating that the meteorological model had more difficulty in reproducing summer transport patterns. Larger fractions of nitrogen compounds were attributed to sources east of the park during the summer than during the spring.

Beem et al. [18], Day et al. [38], and Benedict et al. [3] examined spatial patterns in nitrogen concentrations measured in Northern Colorado during 2006, 2008, and 2009, respectively. Day et al. described patterns in ammonia concentrations in Northeastern Colorado, finding large gradients near sources such as confined animal feeding operations, likely due to dilution and a high dry deposition rate as material moves downwind. Ammonia concentrations in agricultural and urban areas in Colorado east of RMNP were 1-2 orders of magnitude greater than those in the park. The highest ammonia concentrations were in the summer, but mean winter concentrations were higher than spring or fall, likely due to lower mixing heights and wind speeds during the winter. Beem et al. and Benedict et al. found that deposition and airborne concentrations of key reactive nitrogen species were much higher east of the Continental Divide than west and that the highest values measured on the east side of RMNP were associated with easterly flow. Diurnal patterns also showed peak concentrations during the day when easterly upslope flow is most frequent.

Malm et al. [17] described a hybrid deterministic-receptor approach to apportion ammonia measured at RMNP in 2009 during RoMANS II. An inert tracer was simulated by tracking ammonia emissions from 107 source regions with a chemical transport model without deposition or chemistry. Source regions were combined into common transport pathways, and a regression model apportioned the ammonia. On an annual average, about half the ammonia was attributed to sources within Colorado. Of the half from outside the state, the largest fractions were from west rather than east of Colorado, with the largest contributors being California, Utah, and the Snake River valley in Idaho. They hypothesized that diurnal patterns in wind directions and concentrations provide evidence that the mesoscale meteorological model is not adequately capturing the upslope flow at RMNP.

Several themes emerge from this earlier work. The earliest studies focused on oxidized rather than reduced nitrogen, and their respective sources are likely to differ. It is clear that when upslope easterly winds occur on the east side of the Rocky Mountains, they bring emissions from the Front Range and from agricultural areas in Eastern Colorado. So far, despite relatively high ammonia emissions in bordering states to the east, there is no evidence that sources in those states are large contributors to reactive nitrogen in RMNP. The influence of sources to the west is not as well understood. Westerly winds predominate, but easterly winds are more common during precipitation events that lead to wet deposition. Also, the diurnal and seasonal cycle of mountain-valley winds that could bring air pollutants from sources in Central and Eastern Colorado is highly correlated with measured concentrations at RMNP. However, back trajectory analyses indicate that the frequent long-term dribble of emissions from the west may be very important on average, and some studies have hypothesized that the geography of the region is conducive to mixing of emissions from both east and west. The relatively shallow diurnal mountain-valley flow has also proved to be difficult to model, leading to potential overestimation of the influence of sources to the west. Therefore, to accurately determine the influence of sources of nitrogen deposition in RMNP, it is important to be able to accurately model the easterly wind flow when it occurs.

This paper describes the mesoscale meteorological modeling conducted for RoMANS II, including some evaluation of those results and gives details of a back-trajectory-based source apportionment of ammonia concentrations for RoMANS II. The source apportionment technique described in this paper extends the RoMANS I trajectory-based source attribution to include a full year and examines the possibility that virtually moving the modeled receptor site slightly in either space or time might improve the correlation between measured and modeled concentrations by accounting for bias in the modeled wind field.

2. Methodology

2.1. Mesoscale Meteorological Modeling

Gridded meteorological fields for RoMANS II were generated with the weather research and forecasting (WRF) model [39, 40]. The model was run more than dozen times using different physics options, observational meteorological data, observational nudging coefficients, and radii of influence for observational data in an attempt to get the best possible match between model output and observed meteorology. The first runs were with WRF version 3.2, while later runs used version 3.3.1. There were no discernible differences in output attributed to the different versions of the model.

Many model settings were common to all WRF cases. These include the three nested horizontal domains of 36, 12, and 4 km grid sizes identical to those used for the 2006 RoMANS I study [11]. The outer domain, known as the Western Regional Air Partnership (WRAP) domain, covered most of North America. The 12 km domain covered parts of several western states, while the 4 km grid encompassed most of Colorado. North American Regional Reanalysis (NARR) data [41, 42] were used for boundary conditions, initial conditions, and 3-hour analysis nudging on the 36 km domain. Reanalysis rather than higher resolution forecast data were chosen to maximize the observational data in the input. The model top was 100 hPa to match the top of NARR, and there were 34 vertical layers for all runs. Initially, levels were chosen automatically by WRF, but, in later cases, layers were selected manually to increase the resolution near the surface, to match layers used in previous modeling [11] and to correspond to emissions data used for the chemical transport model. Increasing the number of layers near the surface necessitated a factor of 3 reduction in time step (from 180 to 60 seconds) for some days. WRF was run for multiple 3.5-day (84 hours) time periods with overlapping 12-hour spin up times, in part based on the findings of Colle et al. [43], who found that the best precipitation forecasts were 18–36 hours after initialization but that at least 12–18 hours are needed for spin up.

Physics options for initial WRF cases were chosen for simplicity. These are discussed later, but the switch to more carefully considered options did not substantially change the model outcomes. Final choices for the major physics options and reasons for their use are summarized in Table 1.

Table 1: Major physics choices for most WRF model runs, plus reasons. The values in parentheses are the WRF option names in column 1 and the WRF codes for the chosen options in column 2.

Parameterization (WRF namelist variable)
Choice (WRF code)
Reasons

Longwave radiation (ra_lw_physics)
Rapid radiative transfer model (RRTM) (1) [20]
Accurate, fast, accounts for multiple bands, trace gases, and microphysics, commonly used, used in RoMANS I

Shortwave radiation (ra_sw_physics)
Goddard, 2-stream multi-band scheme with ozone from climatology and cloud effects (2) [21, 22]
Recommended in WRF manual for similar scenarios

Surface layer physics (sf_sfclay_physics)
Monin-Obukhov similarity theory provides exchange coefficients to surface scheme. (1) [23, 24]
Compatible with PBL scheme below, WRF manual often recommends, used in RoMANS I

Land surface physics (sf_surface_physics)
Noah land surface model, unified NCEP/NCAR/AFWA scheme with soil temp and moisture in 4 layers, fractional snow cover and frozen soil physics. New features in WRF 3.1 for processes over ice and snow. (2) [25]
Used in NARR and RoMANS I, used operationally at NCEP, compatible with option for PBL scheme below

Planetary boundary layer (PBL) physics (bl_pbl_physics)
Yonsei University scheme, non-local K scheme with explicit entrainment layer and parabolic K profile in unstable mixed layer (1) [26–28]
Used another scheme in RoMANS I that is now an “older scheme”

Findings of Hu et al. (2010) [28] that this had low biases, commonly used

Microphysics (mp_physics)
WRF single-moment 5-class scheme, has mixed-phase and super-cooled water (4) [29]
Relatively cheap, includes super cooled water, snow melt, ice sedimentation, and time-split fall terms, suitable for mesoscale grid sizes

Cumulus clouds (cu_physics)
Kain-Fritsch, deep and shallow convection subgrid scheme with mass flux approach, downdrafts and CAPE removal time scale, not used on 4-km domain (1) [30, 31]
Recommendation of WRF manual, commonly used scheme

The largest differences in modeled output were due to changing the details of the observational nudging on the 4 km domain, including the input data, the nudging coefficients, and the radii of influence. Observational data for nudging on the 4 km domain always included hourly data from a 10-meter tower at the core site and from April 10 to December 27, 2009, from a radar wind profiler in Estes Park, Colorado, approximately 12 km to the NNE of the core site (described in [16] and shown in Figure 1). Early WRF cases also used observational data from the National Center for Atmospheric Research datasets DS461.0 [44] and DS351.0 [45], for surface and upper air data, respectively. Later, data from the Meteorological Assimilation Data Ingest System (MADIS) [46] were used instead because MADIS has a much denser network of observations, particularly along the Front Range. However, as discussed in the following section on model assessment, more data did not always prove to be better. Figure 2 shows the locations of the meteorological observations in the NCAR and MADIS datasets. Panel (b) of Figure 1 also shows the site of a Remote Automated Weather Station (RAWS) in Estes Park. Data from this site are included in the MADIS, but not the NCAR, observational data.

(a)  

(b)

Figure 2: Colorado (450 km north-south by 610 km east-west) with county boundaries and sites of meteorological observations. Panel (a) shows sites in NCAR’s DS461.0 dataset. The red hatched area is RMNP. Site 40 on the eastern edge of RMNP is the RoMANS core site. The remaining sites were grouped as shown for WRF model assessment. Panel (b) shows locations of some MADIS data including locations of the NOAA profiler network (NPN), radiosonde observations (RAOB) and hourly surface (METAR) data.

3. Assessment of WRF Wind Fields

3.1. Means by Modeling Cases: Surface Wind Example

One method to assess WRF’s reproduction of key meteorological variables was to generate many graphs similar to those shown in Figure 3 in which means of modeled and observed values are compared for various spatial and temporal scales. This particularly interesting example shows the east-west component of the wind vector, U, as observed and as modeled from several different cases of WRF. In all cases, WRF was run from November 2008 to November 2009. The top panel shows averages by month for the 45 sites shown in Figure 2(a). The bottom left of Figure 3 has means by hour of day for these 45 sites, and the bottom right shows means by site for all times. In each panel, observed values are in black and modeled values have colored symbols. Results from the first simplest case (red circles) have no observational nudging and mostly default physics options. Blue diamonds show results from the next case, identical to the first, except observational nudging that was added to the 4 km domain using the NCAR data discussed above. Green triangles are from the next attempt, in which the physics options were changed from mostly default values to those shown in Table 1. Magenta triangles are from a case identical to that shown in green, except satellite wind data were removed from the NCAR dataset (discussed below); cyan x-filled squares are from the first case to use MADIS data rather than NCAR for observational nudging, and, finally, the brown squares show results from the final run also using MADIS data but without aircraft observations and with the final observational nudging coefficients.

Figure 3: Mean values of U, the east-west component of the wind vector, for several WRF cases, by month, by hour, and by site. The sites are shown in Figure 2(a). The colored lines in the key represent results from different WRF cases from the earliest and simplest at the end to the final at the beginning. The black line is the observed data.

Examining the data in this way provides insight into potential model bias. From the means by month in the top panel, it is apparent that all WRF cases reproduced the mean seasonal pattern, but for most attempts the mean modeled wind speeds were too high or too westerly for all months. Second, the small differences between the values in blue and those in green on all panels show that changing physics options, on average, had much less impact than changing the input data and other parameters associated with observational nudging. The bottom right panel shows that U is much better modeled at some sites than at others and that site 40, the RoMANS core site, has the poorest reproduction. It is also apparent that while, on average, U is usually over estimated, this is not true for all sites.

The most obvious features in the lower left panel are the odd periodic spikes occurring every 6 hours for two cases. They occur only with observational nudging using NCAR data and are from cases with two different sets of physics options, so it was apparent that the spikes were due to either the observational data or the preprocessing of that data. The only 6-hour data in the NCAR datasets are winds derived from satellites. When these data were removed, the result was the case shown by the magenta triangles. The initial switch from NCAR to MADIS data surprisingly resulted in values shown by the x-filled cyan squares, clearly the worst results. Based on the experience with satellite winds, MADIS aircraft data were eliminated, resulting in the final WRF output used for all further RoMANS II modeling and analyses and shown by the brown squares. There were additional experiments in which the nudging coefficients and radii of influence were varied, but most were for only a single month and for clarity are not shown. Table 2 summarizes the final observational nudging options used in WRF. Other likely options that may have influenced the observational nudging results include the use of the error max and buddy check switches that remove observations when they differ too greatly from the model first-guess field and from neighboring observations. It is possible that in complex terrain too many useful observations are removed by these checks. In this study they were left on, as is the default.

Table 2: Final observational nudging parameters for WRF. Most values are the default values.

Parameter
Value

Domains
4-km only

Variables
Winds, temperature, moisture

Variables nudged in planetary boundary layer
Winds, temperature, moisture

Nudging coefficient for winds
 1/sec (default is )

Nudging coefficient for temperature and moisture
 1/sec

Horizontal radius of influence
60 km (default is 200)

Vertical radius of influence
0.1 eta coordinates

Observational data
10-meter tower at core site, radar wind profiler in Estes Park during Apr–Dec 2009, MADIS (NPN, RAOB, METAR only)

Error max/buddy check
True/True

Temperature: 10/8 K

Horizontal Wind: 13/8 m/sec

RH: 50/40 percent

Sea Level Pressure: 600/800 Pa

Vertical consistency/convective adjustment
False/False

Smoothing
None

Objective analysis scheme
Cressman (Ellipse and Banana variations applied by WRF as warranted). Four passes with default 5, 4, 3, 2 grid cell radii.

Clearly, the WRF results were sensitive to the data used for observational nudging. Initially, it was not clear why aircraft and satellite wind data caused unrealistic increases in wind speeds, but recently Wee et al. [47] showed that bias in the WRF model related to the vertical coordinate system may cause problems when assimilating height-based observations. The graphical analyses do confirm, however, that, in general, observational nudging is beneficial. Mean values of U from the final case are the closest to the observed values on average.

Further spatial and temporal analyses of the modeled winds, for example, contours of wind speed, U and V (the north-south component of the wind vector), and maps of mean observed and modeled wind vectors for different time periods and sites, led to the observation that WRF could correctly generate the diurnal easterly upslope flow over broad areas of Colorado when it occurred but that the modeled easterly winds often did not extend far enough to the west to reach the core site. During time periods when easterly winds were observed there, WRF often had westerly winds at that site but easterly winds in areas farther to the east. Furthermore, as might be expected, WRF’s reproduction of the observed winds was better in areas of less complex terrain, such as the high plains of Eastern Colorado. For example, Figure 4 shows the mean midafternoon wind vectors observed and modeled during June–August 2009. Note the convergence of easterly and westerly winds along the Front Range and the high wind speeds along the Continental Divide. These wind patterns give strength to the hypothesis that this is likely a region where emissions from sources on both sides of the mountains often converge and so have the potential to react.

Figure 4: Mean observed (black) and modeled (red) wind vectors in Colorado for June–August 2009 at 3 : 00 pm local time. Internal lines are county boundaries. The hatched area is Rocky Mountain National Park. Colorado is 450 by 610 km, bounded by 37 and 41 degrees north latitude and by 109 and 102 degrees west longitude.

3.2. Wind Roses of Observed and Final Modeled Winds

Wind roses generated using the Openair Library [48] for R statistical software [49] for the study year for observed and modeled winds from the final WRF run for seven regions in Colorado are shown in Figure 5 (see in Figure 2(a) for a map of the regions). In general, the directional frequency patterns are broadly reproduced by WRF, though the level of accuracy varies by region. As expected, modeled and observed wind directions are most similar in the Northeastern and Southeastern Plains, where the terrain is relatively flat, and perhaps more unexpectedly, also reproduced well in the Southern Mountains. Directional frequencies are most poorly reproduced in the Northern Mountains and the Southwest. The predominant wind direction in the Northern Hemisphere midlatitudes is generally westerly, but, in Colorado, channeling and blocking by terrain make this generality less true in some areas. For example, note the high frequency of north-south winds in the Northern and Southern Front Range, due to channeling along the eastern edge of the Rocky Mountains.

Figure 5: Wind roses for observed and modeled winds in each of seven regions in Colorado from November 2008 to November 2009. Each light gray circle represents 5% frequency. Region locations are shown in Figure 2(a).

Accuracy in the distributions of wind speeds also varies by region. For example, in the Northern Front Range and Northern Mountains, the frequencies of higher-speed winds are underestimated, while frequencies of lower speed winds are overestimated. The reverse is true in the Southwest, while in the Northeast Plains and Southern Mountains the observed wind speed frequencies are more accurate than in other regions.

While, on average, WRF reproduces the winds with reasonable accuracy, at the core site, this was less true. Figure 6 shows wind roses of observed and modeled 10-meter winds at three sites in close proximity to one another: the core site, the radar wind profiler site in Estes Park, and at the RAWS site in Estes Park. These are all for May–November 2009 when the radar wind profiler was operating (see Figure 1 for locations). The observed winds at the core site are channeled by the terrain and arrive most frequently from the northwest with upslope winds mostly from the southeast rather than from due east. In contrast, the modeled output from WRF is much too frequently from the west at this site.

Figure 6: Wind roses showing the percent of winds from each of 12 directional bins for 10-meter winds during May–Nov 2009 as observed and as modeled by WRF at the core site, at the Estes Park radar wind profiler site and at the Estes Park RAWS site. Each light gray circle represents 5% frequency.

At the profiler site, 12 km to the north-northeast, the observed winds are on average rotated about 30 degrees counterclockwise from those at the core site, arriving more frequently from due west and due east, rather than from northwest and southeast. WRF captures this westerly flow fairly accurately but underestimates the frequency of easterly winds. In the adjacent 4 km grid cell, observed winds at a RAWS site also do not have the high frequency of winds from due east that are observed at the profiler. The RAWS site is approximately 3.7 km west of the profiler and about 91 meters higher in elevation. The observed winds at these two locations are quite different, with the RAWS site having predominantly southwesterly rather than westerly flow and with its easterly winds tending to be from the east-northeast rather than the due easterly winds observed at the profiler. As might be expected from adjacent grid cells, the modeled winds for these two locations are similar. Neither of the grid cells has a high enough modeled frequency of easterly winds.

These results led to speculation that, due to the complex terrain, the modeled winds in a nearby grid cell or at a slightly different time might match the true winds at the core site better than those in the spatially correct grid cell. Figure 7 shows wind roses of the modeled winds at the core site and at sites located one 12 km grid cell in each of the cardinal directions. Note that the modeled winds in the grid cell 12 km east are a much better match to the measured winds at the core site (Figure 6) than are the modeled winds in the corresponding WRF grid cell.

Figure 7: Wind roses showing the percent of winds from each of 12 directional bins for 10-meter winds during November 2008 through November 2009 as modeled by WRF at the core site and at sites located 12 km in each of the cardinal directions. Wind speeds are indicated by the colors as shown in the key.

4. Back Trajectories

To further investigate transport of emissions to the core site, hourly ensemble back trajectories with a maximum length of 7 days were generated using version 4.9 revision 346 of the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model [50]. In ensemble mode, HYSPLIT generates 27 trajectories for each start time by using start points at the specified start location, as well as a layer above and a layer below and eight surrounding grid cells. Two sets of HYSPLIT output were generated, one with input from the WRF model discussed above and another using output from the North American mesoscale model at 12 km grid spacing (NAM12) [51, 52] as input. HYSPLIT was started from five locations, the core site and from 12 km north, south, east, and west. Figure 8 shows an example of the 27 ensemble trajectories started from the core site during a single hour. Note that the categorization of an hour into either easterly or westerly flow is not always straightforward. Figure 8(a) shows the horizontal pathways within Colorado, where within the first 1-2 days most ensemble members passed over areas in Northeastern Colorado. However, as Figure 8(b), showing the same trajectories, illustrates, when traced for a longer time, it is evident that these trajectories originated in areas to the west that also have high ammonia emissions, including the Snake River valley in Southern Idaho and agricultural areas in Northern Utah. Figure 8(c) shows the trajectory heights by time. The dashed line is the top of the mixed layer for the central member of the ensemble. Trajectories truncate when reaching the edge of all input data domains. In this example, some ensemble members reached the western edge of the WRF 36 km domain within 4-5 days.

(a) Horizontal trace in Colorado

(b) Horizontal trace in Western United States

(c) Heights by time

Figure 8: Hysplit ensemble back trajectories generated with WRF input for 10 May 2009, at 6 : 00 pm MST. The heaviest line is the trajectory started from the core site. The lighter lines are the remaining 26 ensemble members. In panel (c), the dashed line is the height of the top of the mixed layer for core site trajectory. There is a dot after each day on all panels. Colorado is approximately 610 km across. On panel (a), the distance from the Pacific Coast on the west to the eastern edge of the map is approximately 3200 km.

Figure 9 shows three related analyses generated from the WRF back trajectories started from the core site at 10 m. Statistics in each panel are generated from approximately 256,000 individual trajectories. Figure 9(a) is the overall residence time for the study year. Each endpoint represents an hour of time in a grid cell. It shows where air masses most frequently resided before arriving at the receptor. Westerly flow that tends to bifurcate into either northwesterly or southwesterly predominates. However, there is also easterly transport that can bring emissions from Front Range cities and agricultural areas in eastern Colorado and from states to the east. Figure 9(b) shows the mean concentration of ammonia at the core site when air masses arrived from each grid cell. It is calculated by assigning each trajectory endpoint the value of the measured concentration at the receptor for the trajectory arrival time. The value for each grid cell is the mean of the arrival concentrations for endpoints associated with trajectories that passed through that cell. This panel illustrates that when air masses arrive from the east the concentration at the receptor is higher than when air masses arrive from the west. However, it does not take into account transport frequency. By combining the information contained in both top panels, it can be seen that transport from the west is common, but air masses from this direction tend to arrive with relatively lower ammonia concentrations. The inverse is true for air masses from the east; they are less frequent but on average arrive with higher concentrations. These statistics are combined to generate the map in Figure 9(c), which is generated by multiplying the results of the first two panels to create concentration-weighted residence times. These results take into account the mean concentration as well as the transport frequency and can be thought of as an indication of the cumulative impact of a region. The pattern is similar to the residence time panel, with some subtle differences. Areas to the east are weighted relatively higher in Figure 9(c) than in Figure 9(a). For example, the high average concentrations during times when trajectories arrived from Eastern Texas cause values in that region to be higher in the concentration-weighted residence time than in the simple residence time. The color break points were chosen to best illustrate the regional differences in each statistic. Note from the geometric scales on Figures 9(a) and 9(c) compared to the linear scale on Figure 9(b) that differences in transport frequency are much larger than differences in mean concentration and so dominate the results in Figure 9(c). These results indicate that sources to both the east and west are important contributors to ammonia in RMNP.

(a) Residence time

(b) Mean concentration at receptor

(c) Concentration-weighted residence time

Figure 9: Three analyses of WRF back trajectories for 24-hour measured ammonia at the core site for the study year shown on maps of the Western United States with state boundaries. (a) Overall residence time or likelihood of an air mass arriving from each grid cell. (b) Mean ammonia concentration at the receptor when air mass arrived from each grid cell. (c) Concentration-weighted residence time. The distance from the Pacific Coast to the eastern edge of each map is approximately 2500 km.

5. Measured Concentrations

The measured concentrations obtained as part of RoMANS II are described in detail in Malm et al. [17] and Benedict et al. [2, 3]. The 24-hour data from the URG denuder are considered the most reliable ammonia measurement and are the data used for TrMB modeling. Figure 10 shows a timeline of these concentrations observed at the core site. In general, concentrations are higher during the warmer months of the year. Values above the 90th percentile occurred during all months from March to September.

Figure 10: Observed 24-hour average concentrations of ammonia at the core site from November 2008 to November 2009. The horizontal lines are the 10th, 50th, and 90th percentiles. Many peak values are labeled with month/day.

6. Trajectory Mass Balance Model

The trajectory mass balance model (TrMB) has been previously described in detail (e.g., [16, 53, 54]). It is a linear-regression-based receptor model with measured concentrations at the receptor as the dependent variable and counts of back trajectory endpoints in selected source regions as independent variables. Trajectory endpoint counts are interpreted as time spent in the source region prior to arrival at the receptor. Regression coefficients account for mean emissions, dispersion, chemistry, and deposition. The coefficient multiplied by the endpoints is an estimate of the mean contribution from a source region over the time period of the observations.

TrMB was used to estimate the mean seasonal contributions from each of the source regions shown in Figure 11 to the measured ammonia concentrations at the core site. In addition, tests of the sensitivity to the timing and spatial accuracy of the winds were conducted in which trajectories from the preceding and following day and from trajectory start points in grid cells 12 km east, north, west, and south were compared to results from using trajectories started in the correct grid cell and on the correct day.

Figure 11: Source areas used for TrMB modeling and ammonia emissions [32] used for RoMANS II chemical transport modeling on the 36 km (a) and 4 km domains (b). The distance from the eastern edge of the map to the Pacific Coast in panel (a) is approximately 4000 km. In panel (b), Colorado is approximately 610 km from east to west.

For collinear source regions, the endpoint counts were summed to create one combined, rather than two individual, regions. The test was iterative so that combined source regions were retested against the remaining regions until no collinear areas remained. The attributions of the regions were then reapportioned back to the original smaller regions using a weighting factor of mean emissions/mean distance. The definition of collinear is arbitrary; for this study, trials using 70%, 80%, and 90% correlation were used as the limits for collinearity. The number of collinear source regions varied depending on this cutoff value, as well as the distributions of the endpoints between the regions, which depended on season, input meteorology, start location, and so forth. Of the 35 source regions, the median numbers of collinear regions for cutoff values of 70%, 80%, and 90% are 13 (range 2–23), 7 (0–17), and 2 (0–12), respectively.

For each of the four seasons, winter (December, January, and February), spring (March, April, and May), summer (June, July, and August), and fall (September, October, and November), 1080 sensitivity trials of the TrMB model were run by varying four maximum trajectory heights (1000 meters, 3000 meters, top of the mixed layer, and no limit), three trajectory lengths (3, 5, and 7 days), three concentration times (as reported, one day early, and one day late), two meteorological input data sets (WRF and NAM12), five start locations (core site, 12 km east, west, north, and south), and three collinearity limits (70%, 80%, and 90%).

7. Results

One measure of model performance is to examine how well TrMB is able to match the temporal characteristics in the observed data. A better match, as signified by a higher coefficient of variation between modeled and measured concentrations, does not ensure an accurate model, but a poorer match gives lower confidence in the results. The median values for varying subsets of TrMB trials are shown in Table 3. The two biggest determinants of for the RoMANS II ammonia concentrations are season and the cutoff value for collinearity. Seasonally, the highest median (0.54) is for fall and the lowest (0.39) is for summer. Spring (0.40) was only slightly better than summer, and winter (0.49) was almost as good as fall. The higher the correlation cutoff for collinearity is, the higher the median value is. This is because the more the source regions are combined, the lower the degrees of freedom in the regression are. There is a small effect on from moving the receptor 12 km in any direction. The best fit was for 12 km east (0.52) and the worst was for 12 km west (0.44). This is another indication that the easterly winds at the core site may be underestimated. Concentrations on the correct day or a day early had similar median , but moving the concentrations to a day late improved the performance for unknown reasons. Trajectory duration had little effect on , indicating that most source areas are usually less than 5 days transport time from RMNP. In general, the higher the maximum trajectory endpoint height is, the higher the , but this effect was also small. Contrary to what might be expected due to the difference in grid size, use of NAM12 rather than WRF as input to HYSPLIT resulted in a better fit. This is consistent with results from RoMANS I [16].

Table 3: Median values of , percent of ammonia attributed to Colorado and percent attributed to Northeast Colorado for all TrMB trials. Total number of trials is 4320, so trials/row is 4320 divided by the number of categories. For example there are four seasons, so trials/season is .

Category

% from Colorado
% from NE Colorado

By season
Winter
0.49
47
5

Spring
0.40
44
9

Summer
0.39
52
16

Fall
0.54
46
17

By start location
Core site
0.45
48
13

12 km east
0.52
48
17

12 km west
0.44
48
8

12 km north
0.47
47
12

12 km south
0.49
47
11

By maximum trajectory length
3 days
0.42
52
12

5 days
0.45
45
12

7 days
0.44
44
12

By maximum endpoint height allowed
1 km
0.41
52
13

Top of MH
0.43
52
13

3 km
0.44
44
11

No limit
0.46
42
11

By input meteorology
WRF
0.41
46
10

NAM12
0.46
48
14

By concentration timing
Correct time
0.42
48
14

1 day early
0.40
46
10

1 day late
0.48
49
11

By R cutoff for collinearity
0.7
0.36
53
11

0.8
0.44
48
13

0.9
0.50
43
12

Table 3 also shows median source attribution results for all of Colorado and for Northeast Colorado (sources 6, 7, and 11 in Figure 11) for variations in input. These regions were chosen to facilitate examining sensitivity of apportionment to source areas that were close versus far and to the east versus west. The median attributions for the sum of all source regions within the state are quite stable. No matter which parameter is varied, it remains in the 42–53% range. By season, the greatest attribution to Colorado was 52% in the summer, with the remaining seasons all having median attributions of 44–47%. As expected, restricting the trajectories to shorter durations or lower heights somewhat increases the attributions to in-state sources because these changes cause relatively fewer endpoints to reside in distant source regions. Reducing the correlation cutoff for collinearity also increases the attribution to in-state sources. This is due to the distance weighting used when splitting the combined sources apart for the final attribution.

The final sensitivity test, used to test the impact of easterly flow, was attributions to Northeast Colorado. Median values range from 5% to 17%. The largest differences were, by season, ranging from 5% in the winter to 17% in the fall. Attributions for summer (median 16%) were more similar to fall, and spring (9%) was more similar to winter. Moving the start point 12 km west gave a lower median attribution at 8% than moving it east, which resulted in a 17% attribution. Restricting the trajectories in length made no difference in attribution to Northeast Colorado, while restricting the height gives slightly higher attributions to this region because the number of endpoints in more distant regions is reduced, but the differences are very small.

Table 4 and Figure 12 give the final TrMB source attribution results using concentrations only at the correct time but allowing all other input parameters discussed above to vary. Table 4 has seasonal results for aggregated regions and for comparison also shows results from the 2006 RoMANS I study [16], which included only spring and summer measurement campaigns. See Table 5 and Figure 11 for the source areas in each region. Annual concentration-weighted means of the seasonal attributions are also shown in the bottom row of Table 4 and, for the individual source areas, in the numerators of the numbers in Figure 12. The denominators in Figure 12 are the annual mean percent of back trajectory endpoints in each area rounded to the nearest whole percent. It is interesting to compare the source attributions to the percent of endpoints for the various regions. For most source areas within Colorado, the attribution is the same or greater than the percent of endpoints. The exceptions are two areas in Northwestern Colorado from which there is frequent transport, but emissions are relatively low. TrMB results are thus consistent with the expectation that, due to deposition, chemical reaction, and dispersion, closer sources contribute more per endpoint than do more distant regions. Also, source area 6, the area just to the east of the receptor, contributes much more per endpoint than any other regions. This area has high emissions of ammonia due to confined animal feeding operations and other agricultural activities. The TrMB model does not explicitly include emissions data as input (except as a weighting factor in the case of collinear regions), yet the results do reflect the higher emissions in that region. For sources outside of Colorado, with a few exceptions, source areas to the east of the state have a higher percent attribution of ammonia than percent of trajectory endpoints, while the reverse is usually true for areas to the west. This again reflects the higher emissions in areas to the east. The higher mean attributions for sources to the west are due to the much greater frequency of transport from those regions.

Table 4: Mean measured 24-hour ammonia concentrations by season in the far right column and percent attributions ± standard deviations by TrMB analysis for the 2006 RoMANS I (italic font) and 2009 RoMANS II studies. The values in the bottom row (bold font) are concentration weighted means of the four 2008/09 seasons.The RoMANS II results are composite results from 360 TrMB trials for each season using concentrations for the correct time, but allowing all other parameters shown in Table 3 to vary. See Table 5 and Figure 11 for locations.

Season
NE Colorado (%)
SE Colorado (%)
CO Front Range (%)
Western Colorado (%)
Total Colorado (%)
Eastern U S

(%)
Western U S

(%)
URG 24-hr Conc (g/m³)

RoMANS I

spring 2006
7 ± 1
NA
9 ± 3
39 ± 4
56 ± 5
8 ± 2
35 ± 4
0.14

RoMANS I

summer 2006
9 ± 5
NA
15 ± 4
24 ± 5
50 ± 7
21 ± 3
29 ± 5
0.35

RoMANS II

winter 2008/09

0.09

RoMANS II

spring 2009

0.32

RoMANS II

summer 2009

0.31

RoMANS II

fall 2008/09

0.23

RoMANS II

overall 2008/09
15 ± 8
2 ± 4
5 ± 4
26 ± 11
48 ± 12
9 ± 8
43 ± 11
0.24

Table 5: Individual source areas included in each of the regions in Table 4. See Figure 11 for a map of the regions.

Region
Source areas in Figure 11

Northeast Colorado
6, 7, 11

Southeast Colorado
10

Colorado Front Range
3, 5

Western Colorado
1, 2, 4, 8, 9, 12

Eastern Colorado
3, 5, 6, 7, 10, 11

All Colorado
1–12

Eastern U S
22–27, 34, 35

Western U S
13–21, 28–33

Figure 12: Source areas for TrMB modeling of ammonia. Top left is the Western US; the inset on the lower right is Colorado. County boundaries and the receptor site are shown in the inset. The numbers are the annual mean TrMB percent ammonia attribution in the numerator calculated as in the bottom row of Table 4 and the percent of back trajectory endpoints in the denominator. Values are rounded to the nearest integer; all areas had at least a fraction of a percent of the total endpoints.

TrMB estimates that approximately half of the ammonia at the core site is from sources within Colorado, though winter and summer contribute higher percentages than spring and fall. Areas to the west contribute more on average than areas to the east. Western Colorado contributes more than Eastern Colorado except during the fall, and Western United States contributes much more than Eastern United States during all seasons. Of areas within Colorado, Western Colorado contributes 18–34% on average, depending on season; Northeast Colorado adds 10–17%, the Front Range adds 2–8%, and Southeast Colorado adds 0–3%. For areas outside the state, Western United States contributes 39–47%, while states to the east add only 5–12%. The 2008-2009 results vary from those for spring and summer 2006, but there are a few consistent patterns. Sources in Western United States and Northeast Colorado contributed more in both seasons of 2009, while sources in Eastern United States and the Colorado Front Range contributed more in both seasons of 2006. The measured mean concentrations were similar during the two summers, but the mean spring concentration was higher by about a factor of 2 during 2009 as in 2006.

8. Discussion and Conclusions

Previous studies of sources of nitrogen in Rocky Mountain National Park and surrounding areas have shown that concentrations and deposition are generally higher when there are upslope easterly winds and that easterly winds are more likely during precipitation events (wet deposition) than during dry weather. There are also indications from this study and previous work that this is a convergence zone where emissions from both east and west of the Continental Divide meet and thus can potentially react. Also, simple maps of back trajectories show that categorizing times into purely easterly or westerly flow is too simplistic because curvature in the transport pathways can simultaneously bring emissions from both directions.

The complex terrain in the mountains makes meteorological modeling more difficult than in areas of flatter terrain. Small-scale valley and slope winds are not adequately modeled at the 4 km grid scale within and near RMNP where mean observed wind directions vary substantially within a few kilometers. Observational nudging improved the WRF model performance, but only when satellite winds and aircraft observations were removed from the observational data. WRF-modeled winds were more sensitive to observational nudging options than to physics options.

Statistical analyses of back trajectories for the 2008-2009 RoMANS II study corroborate earlier studies that show that transport from the west predominates, but when there is flow from the eastern half of the state, ammonia concentrations are higher than that when the transport pathway is more purely from the west. Despite high ammonia emissions in states east of Colorado, simulated transport from these states is infrequent, and on average they were found to have small contributions to the measured ammonia concentrations at RMNP.

Use of the TrMB model to estimate source attributions for 24-hour ammonia concentrations measured at RMNP was completed for the RoMANS II study. Experiments conducted to determine if the model was sensitive to moving the receptor location slightly were prompted by the observation that easterly winds were not generated by the WRF model as frequently as they were observed at the core site and that modeled winds 12 km east of the core site matched the observations better than modeled winds at the actual site. The TrMB source attribution results were often not very sensitive to moving the start location or start time, especially for areas outside of Colorado. Moving the trajectory start point to the east results in attributions to source regions in Eastern Colorado that are a few percentage points higher, and the ability of the model to reproduce the temporal patterns in the measured concentrations is also slightly better than when the trajectory is started from the actual measurement site. TrMB results show that about half the ammonia measured at the core site is from sources within the state and that, on average, sources to the west contribute more than sources to the east. This was true for all seasons. Ammonia deposits rapidly and is readily absorbed by cloud and rain drops, so the typical atmospheric lifetime is several hours to a few days [55]. If sources as distant as Idaho and California contribute to observed gaseous ammonia concentrations at RMNP, either the ammonia was lofted high enough to have little interaction with the surface during transport or it converted to a related compound such as ammonium nitrate, which has a much slower deposition velocity and then back to ammonia when close to the measurement location.

TrMB is only one piece of several source attribution techniques employed for RoMANS II.

Malm et al. [17] reported the results from a hybrid deterministic-receptor technique. On an annual average, this technique also apportioned about half the ammonia to sources within Colorado. Of the half from outside the state, the largest fractions were from west rather than east of Colorado. However, comparison with TrMB results is difficult due to different source region definitions.

Results from chemical transport modeling have not yet been published and reconciliation of the various source attribution methodologies is ongoing.

Disclaimer

The assumptions, findings, conclusions, judgments, and views presented herein are those of the authors and should not be interpreted as necessarily representing official National Park Service policies.

Conflict of Interests

The authors declare that they have no conflict of interests regarding the publication of this paper.

References

W. M. Lewis Jr., M. C. Grant, and J. F. Saunders III, “Chemical patterns of bulk atmospheric deposition in the state of Colorado,” Water Resources Research, vol. 20, no. 11, pp. 1691–1704, 1984.

K. B. Benedict, S. M. Kreidenweis, B. A. Schichtel, W. C. Malm, and J. L. Collett Jr., “A seasonal nitrogen deposition budget for Rocky Mountain National Park,” Ecological Applications, vol. 23, no. 5, pp. 1156–1169, 2013.

K. B. Benedict, D. E. Day, F. M. Schwandner et al., “Observations of atmospheric reactive nitrogen species in Rocky Mountain National Park and across Colorado,” Atmospheric Environment, vol. 64, pp. 66–76, 2013.

D. A. Burns, “Atmospheric nitrogen deposition in the Rocky Mountains of Colorado and southern Wyoming—a review and new analysis of past study results,” Atmospheric Environment, vol. 37, no. 7, pp. 921–932, 2003.

C. M. B. Lehmann, V. C. Bowersox, and S. M. Larson, “Spatial and temporal trends of precipitation chemistry in the United States, 1985–2002,” Environmental Pollution, vol. 135, no. 3, pp. 347–361, 2005.

D. H. Campbell, C. Kendall, C. C. Y. Chang, S. R. Silva, and K. A. Tonnessen, “Pathways for nitrate release from an alpine watershed: Determination using delta N-15 and delto O-18,” Water Resources Research, vol. 38, no. 5, pp. 10-1–10-9, 2002.

J. S. Baron, T. M. Schmidt, and M. D. Hartman, “Climate-induced changes in high elevation stream nitrate dynamics,” Global Change Biology, vol. 15, no. 7, pp. 1777–1789, 2009.

J. S. Baron, H. M. Rueth, A. M. Wolfe et al., “Ecosystem responses to nitrogen deposition in the Colorado Front Range,” Ecosystems, vol. 3, no. 4, pp. 352–368, 2000.

M. E. Fenn, J. S. Baron, E. B. Allen et al., “Ecological effects of nitrogen deposition in the western United States,” BioScience, vol. 53, no. 4, pp. 404–420, 2003.

M. E. Fenn, R. Haeuber, G. S. Tonnesen et al., “Nitrogen emissions, deposition, and monitoring in the western United States,” BioScience, vol. 53, no. 4, pp. 391–403, 2003.

M. A. Rodriguez, M. G. Barna, K. A. Gebhart et al., “Modeling the fate of atmospheric reduced nitrogen during the Rocky mountain atmospheric nitrogen and sulfur study (RoMANS): performance evaluation and diagnosis using integrated process analysis,” Atmospheric Environment, vol. 45, no. 1, pp. 223–234, 2011.

D. D. Parrish, “Systematic variations in the concentration of NOx (NO plus NO₂) at Niwot Ridge, Colorado,” Journal of Geophysical Research, vol. 95, no. 2, pp. 1817–1836, 1990.

J. Baron and A. S. Denning, “The influence of mountain meteorology on precipitation chemistry at low and high elevations of the Colorado Front Range, U.S.A,” Atmospheric Environment A General Topics, vol. 27, no. 15, pp. 2337–2349, 1993.

M. Losleben, N. Pepin, and S. Pedrick, “Relationships of precipitation chemistry, atmospheric circulation, and elevation at two sites on the Colorado front range,” Atmospheric Environment, vol. 34, no. 11, pp. 1723–1737, 2000.

W. C. Malm, J. L. Collett Jr., M. G. Barna et al., Rocky Mountain Atmospheric Nitrogen and Sulfur Study (RoMANS) Final Report, 2009, http://nature.nps.gov/air/Pubs/regionPark.cfm.

K. A. Gebhart, B. A. Schichtel, W. C. Malm, M. G. Barna, M. A. Rodriguez, and J. L. Collett, “Back-trajectory-based source apportionment of airborne sulfur and nitrogen concentrations at Rocky Mountain National Park, Colorado, USA,” Atmospheric Environment, vol. 45, no. 3, pp. 621–633, 2011.

W. C. Malm, B. A. Schichtel, M. G. Barna et al., “Aerosol species concentrations and source apportionment of ammonia at Rocky Mountain National Park,” Journal of the Air & Waste Management Association, vol. 63, no. 11, pp. 1245–1263, 2013.

K. B. Beem, S. Raja, F. M. Schwandner et al., “Deposition of reactive nitrogen during the Rocky Mountain Airborne Nitrogen and Sulfur (RoMANS) study,” Environmental Pollution, vol. 158, no. 3, pp. 862–872, 2010.

N. J. Doesken, R. A. Pielke Sr., and O. A. P. Bliss, Climate of Colorado, Climatography of the United States No. 60, Colorado Climate Center, Atmospheric Science Department, Colorado State University, Fort Collins, Colorado, 2003.

E. J. Mlawer, S. J. Taubman, P. D. Brown, M. J. Iacono, and S. A. Clough, “Radiative transfer for inhomogeneous atmospheres: RRTM, a validated correlated-k model for the longwave,” Journal of Geophysical Research D: Atmospheres, vol. 102, no. 14, pp. 16663–16682, 1997.

M. D. Chou, “Parameterizations for the absorption of solar radiation by O₂ and CO₂ with application to climate studies,” Journal of Climate, vol. 3, pp. 209–217, 1990.

M. D. Chou, “A solar radiation model for use in climate studies,” Journal of the Atmospheric Sciences, vol. 49, no. 9, pp. 762–772, 1992.

A. S. Monin and A. M. Obukhov, “Osnovnye zakonomernosti turbulentnogo peremeshivanija v prizemnom sloe atmosfery (Basic laws of turbulent mixing in the atmosphere near the ground),” Trudy Geofizicheskogo Instituta, Akademiya Nauk SSSR, vol. 24, no. 151, pp. 163–187, 1954.

T. Foken, “50 years of the Monin-Obukhov similarity theory,” Boundary-Layer Meteorology, vol. 119, no. 3, pp. 431–447, 2006.

M. B. Ek, K. E. Mitchell, Y. Lin et al., “Implementation of Noah land surface model advances in the National Centers for Environmental Prediction operational mesoscale Eta model,” Journal of Geophysical Research D: Atmospheres, vol. 108, no. 22, pp. 1–16, 2003.

S.-Y. Hong, Y. Noh, and J. Dudhia, “A new vertical diffusion package with an explicit treatment of entrainment processes,” Monthly Weather Review, vol. 134, no. 9, pp. 2318–2341, 2006.

S.-Y. Hong and S.-W. Kim, “Stable boundary layer mixing in a vertical diffusion scheme,” in Proceedings of the 9th Annual WRF User's Workshop, National Center for Atmospheric Research, 2008.

X.-M. Hu, J. W. Nielsen-Gammon, and F. Zhang, “Evaluation of three planetary boundary layer schemes in the WRF model,” Journal of Applied Meteorology and Climatology, vol. 49, no. 9, pp. 1831–1844, 2010.

S.-Y. Hong, J. Dudhia, and S.-H. Chen, “A revised approach to ice microphysical processes for the bulk parameterization of clouds and precipitation,” Monthly Weather Review, vol. 132, no. 1, pp. 103–120, 2004.

J. S. Kain and J. M. Fritsch, “Convective parameterization for mesoscale models: the Kain-Fritch scheme,” in The Representation of Cumulus Convection in Numerical Models, K. A. Emanuel and D. J. Raymond, Eds., American Meteorological Society, 1993.

J. S. Kain and J. Kain, “The Kain—Fritsch convective parameterization: an update,” Journal of Applied Meteorology, vol. 43, no. 1, pp. 170–181, 2004.

Z. Adelman and M. Omary, Emissions Modeling Final Report—Developing 2009 Emissions For the NPS-ARD RoMANS Study, Chapel Hill: Institute for the Environment, University of North Carolina.

W. D. Neff, “The Denver brown cloud studies from the perspective of model assessment needs and the role of meteorology,” Journal of the Air and Waste Management Association, vol. 47, no. 3, pp. 269–285, 1997.

H. Sievering, D. Rusch, and L. Marquez, “Nitric acid, particulate nitrate and ammonium in the continental free troposphere: nitrogen deposition to an alpine tundra ecosystem,” Atmospheric Environment, vol. 30, no. 14, pp. 2527–2537, 1996.

K. Heuer, K. A. Tonnessen, and G. P. Ingersoll, “Comparison of precipitation chemistry in the Central Rocky Mountains, Colorado, USA,” Atmospheric Environment, vol. 34, no. 11, pp. 1713–1722, 2000.

G. P. Ingersoll, K. A. Tonnessen, D. H. Campbell, B. R. Glass, and A. O. Torizzo, “Effect of storm trajectories on snowfall chemistry in Rocky Mountain National Park, Colorado,” in Proceedings of the 69th Annual Meeting Western Snow Conference, pp. 32–42, April 2001.

H. Sievering, T. Kelly, G. McConville, C. Seibold, and A. Turnipseed, “Nitric acid dry deposition to conifer forests: niwot Ridge spruce-fir-pine study,” Atmospheric Environment, vol. 35, no. 22, pp. 3851–3859, 2001.

D. E. Day, X. Chen, K. A. Gebhart et al., “Spatial and temporal variability of ammonia and other inorganic aerosol species,” Atmospheric Environment, vol. 61, pp. 490–498, 2012.

The Weather Research and Forecasting Model, 2012, http://wrf-model.org/index.php

W. C. Skamarock, J. B. Klemp, J. Dudhia et al., A Description of the Advanced Research WRF Version 3, 2008, NCAR Technical Note NCAR/TN-475+STR.

F. Mesinger, G. DiMego, E. Kalnay et al., “North American regional reanalysis,” Bulletin of the American Meteorological Society, vol. 87, no. 3, pp. 343–360, 2006.

North American Regional Reanalysis home page, 2009, http://wwwt.emc.ncep.noaa.gov/mmb/rreanl/index.html

B. A. Colle, C. F. Mass, and K. J. Westrick, “MM5 precipitation verification over the Pacific Northwest during the 1997-99 cool seasons,” Weather and Forecasting, vol. 15, no. 6, pp. 730–744, 2000.

UCAR, NCEP ADP Global Surface Observational Weather Data, October 1999, continuing, http://rda.ucar.edu/datasets/ds461

UCAR, NCEP ADP Global Upper Air Observational Weather Data, October 1999, continuing, http://rda.ucar.edu/datasets/ds351.0/

“National Oceanic and Atmospheric Administration (NOAA_,” Meteorological Assimilation Data Ingest System (MADIS), 2010, http://madis.noaa.gov/

T.-K. Wee, Y.-H. Kuo, D.-K. Lee, Z. Liu, and W. Wang, “Two overlooked biases of the Advanced Research WRF (ARW) model in geopotential height and temperature,” Monthly Weather Review, vol. 140, no. 12, pp. 3907–3918, 2012.

D. C. Carslaw and K. Ropkins, “Openair—an r package for air quality data analysis,” Environmental Modelling and Software, vol. 27-28, pp. 52–61, 2012.

R Development Core Team, A Language and Environment For Statistical Computing, R Foundation For Statistical Computing, Vienna, Austria, 2011, http://www.R-project.org.

R. R. Draxler and G. D. Hess, “An overview of the HYSPLIT_4 modelling system for trajectories, dispersion and deposition,” Australian Meteorological Magazine, vol. 47, no. 4, pp. 295–308, 1998.

Z. I. Janjic, “A nonhydrostatic model based on a new approach,” Meteorology and Atmospheric Physics, vol. 82, no. 1–4, pp. 271–285, 2003.

National Weather Service Environmental Modeling Center, North American Mesoscale Forecast System (NAM), 2012, http://www.emc.ncep.noaa.gov/NAM/.php.

M. Pitchford and A. Pitchford, “Analysis of regional visibility in the southwest using principal component and back trajectory techniques,” Atmospheric Environment A General Topics, vol. 19, no. 8, pp. 1301–1316, 1985.

K. A. Gebhart, B. A. Schichtel, M. G. Barna, and W. C. Malm, “Quantitative back-trajectory apportionment of sources of particulate sulfate at Big Bend National Park, TX,” Atmospheric Environment, vol. 40, no. 16, pp. 2823–2834, 2006.

R. D. Saylor, E. S. Edgerton, B. E. Hartsell, K. Baumann, and D. A. Hansen, “Continuous gaseous and total ammonia measurements from the southeastern aerosol research and characterization (SEARCH) study,” Atmospheric Environment, vol. 44, no. 38, pp. 4994–5004, 2010.