Assessing the roles emission sources and atmospheric processes play in simulating δ¹⁵N of atmospheric NO_x and NO₃⁻ using CMAQ (version 5.2.1) and SMOKE (version 4.6)

Nitrogen oxides (NO_x= nitric oxide (NO) + nitrogen dioxide (NO₂)) are important trace gases that affect atmospheric chemistry, air quality, and climate. Contemporary development of NO_x emissions inventories is limited by the understanding of the roles of vegetation (net NO_x source or net sink), vehicle emissions from gasoline- and diesel-powered vehicles, the application of NO_x emission control technologies, and accurate verification techniques. The nitrogen stable isotope composition (δ¹⁵N) of NO_x is an effective tool to evaluate the accuracy of the NO_x emission inventories, which are based on different assumptions. In this study, we traced the changes in δ¹⁵N values of NO_x along the “journey” of atmospheric NO_x, driven by atmospheric processes after different sources emit NO_x into the atmosphere. The ¹⁵N was incorporated into the emission input dataset, generated from the US EPA trace gas emission model SMOKE (Sparse Matrix Operator Kernel Emissions). Then the ¹⁵N-incorporated emission input dataset was used to run the CMAQ (Community Multiscale Air Quality) modeling system. By enhancing NO_x deposition, we simulated the expected δ¹⁵N of NO${}_{\mathrm{3}}^{-}$, assuming no isotope fractionation during chemical conversion or deposition. The simulated spatiotemporal patterns in NO_x isotopic composition for both SMOKE outputs (simulations under the “emission only” scenario) and CMAQ outputs (simulations under the “emission + transport + enhanced NO_x loss” scenario) were compared with corresponding measurements in West Lafayette, Indiana, USA. The simulations under the emission + transport + enhanced NO_x loss scenario were also compared to δ¹⁵N of NO${}_{\mathrm{3}}^{-}$ at NADP (National Atmospheric Deposition Program) sites. The results indicate the potential underestimation of emissions from soil, livestock waste, off-road vehicles, and natural-gas power plants and the potential overestimation of emissions from on-road vehicles and coal-fired power plants, if only considering the difference in NO_x isotopic composition for different emission sources. After considering the mixing, dispersion, transport, and deposition of NO_x emission from different sources, the estimation of atmospheric δ¹⁵N(NO_x) shows better agreement (by ∼ 3 ‰) with observations.

1 Introduction

Nitrogen oxides (NO_x= NO + NO₂) are important trace gases that affect atmospheric chemistry, air quality, and climate. The main sources of tropospheric NO_x are anthropogenic emissions from vehicles, power plants, agriculture, livestock waste, as well as natural emissions from lightning and the by-product of nitrification and denitrification occurring in soil (Galloway et al., 2004). The NO_x photochemical cycle generates OH and HO₂ radicals, organic peroxy radicals (RO₂), and ozone (O₃), which ultimately oxidize NO_x into NO_y (NO_y = NO_x+ HONO + HNO₃+ HNO₄+ N₂O₅+ other N oxides). During the photochemical processes that convert NO_x to NO_y, ground-level concentrations of O₃ become elevated, and secondary particles are generated (Seinfeld and Pandis, 2006). Secondary aerosols are hazardous to human health (Lighty et al., 2000) and affect cloud physics, enhancing the reflection of solar radiation (Schwartz, 1996). Thus, the importance of NO_x in air quality, climate, and human and environmental health makes understanding the spatial and temporal variation in the sources of NO_x a vital scientific question.

Despite years of research, however, there are still several significant uncertainties in the NO_x budget. About 15 %–40 % of global NO_x emissions, ranging from 4 to 15 Tg N yr⁻¹, is derived from global soil NO_x emissions, yet evaluating and verifying emission rates using laboratory measurements, field measurements, and satellite observations is still a challenge (Jaegleì et al., 2005; Yan et al., 2005; Stehfest and Bouwman, 2006; Vinken et al., 2014; Rasool et al., 2016). Soil NO_x emissions vary by different biome types, meteorological conditions, N fertilizer application, and soil physicochemical properties (Ludwig et al., 2001). Furthermore, the role of vegetation is to act as a net source of atmospheric NO_x when ambient NO_x concentration is below the “compensation point” versus acting as a net sink of atmospheric NO_x when ambient NO_x concentrations are above it (Johansson, 1987; Thoene et al., 1996; Slovik et al., 1996; Weber and Rennenberg, 1996). This significantly impacts the biotic NO_x emission inventory (Almaraz et al., 2018). Uncertainties also exist in the amount of NO_x emitted during the combustion of fossil fuels by vehicles and industry. According to Parrish (2006), the estimation of on-road vehicle NO_x emission has at least 10 % to 15 % uncertainty. For the mileage-based algorithm, which is used in the National Emission Inventory (NEI), the uncertainty is caused by the limited number of sites to determine the emission factors of vehicle classifications and emission types (Ingalls, 1989; Pierson et al., 1990; Fujita et al., 1992; Pierson et al., 1996; Singer and Harley, 1996). The uncertainty in power plant NO_x emissions results from the choice of emission control technologies, of which the removal efficiencies of NO_x emission are different. NO_x removal by low NO_x burning, over-fire air reduction, and selective non-catalytic reduction is highly variable, ranging from 50 % to 75 % (Srivastava et al., 2005).

The nitrogen stable isotope composition (δ¹⁵N) of NO_x might be a useful tool to help resolve the uncertainties of how NO_x emission sources vary in space and time because NO_x sources have distinctive ¹⁵N $/$ ¹⁴N ratios (Ammann et al., 1999; Felix et al., 2012; Felix and Elliott, 2013; Fibiger et al., 2014; Heaton, 1987; Hoering, 1957; Miller et al., 2017; Walters et al., 2015a, b, 2018). This variability in NO_x ¹⁵N $/$ ¹⁴N ratios is quantified by

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathit{\delta }}^{\mathrm{15}}{\mathrm{NO}}_x\left(\mathrm{‰}\right)& =\left[\left({}^{\mathrm{15}}{\mathrm{NO}}_x{/}^{\mathrm{14}}{\mathrm{NO}}_x\right)/{\left({}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{2}}{/}^{\mathrm{14}}{\mathrm{N}}_{\mathrm{2}}\right)}_{\mathrm{air}}-\mathrm{1}\right]\\ & \times \mathrm{1000})\end{array},\end{array}$$

where ¹⁵NO_x $/$ ¹⁴NO_x is the measurement of ¹⁵N $/$ ¹⁴N in atmospheric NO_x, compared with the ratios in air N₂=0.0036 (for brevity, the δ¹⁵N value of any NO_y compound will be denoted as δ¹⁵NO_y, e.g., δ¹⁵NO${}_{\mathrm{3}}^{-}$). Previous research has shown that there are unique differences in δ¹⁵N values for NO_x from different emission sources and significant variations within each source (Fig. 1). This uniqueness can potentially be used to partition the relative importance of various NO_x sources in a mixed atmosphere. For example, Redling et al. (2013) found higher δ¹⁵N of NO₂ in samples collected closer to the highway compared to those adjacent to a forest, showing the emissions from vehicles were dominant near the highway. A strong positive correlation between δ¹⁵NO${}_{\mathrm{3}}^{-}$ and NO_x emission from coal-fired power plants within 400 km radial area of study sites of deposition suggests local power plant NO_x emissions impacted regional NO_x budgets (Elliott et al., 2007, 2009). What is lacking is a systematic way of evaluating δ¹⁵NO_y values in numerous studies in the context of NO_x sources, regional emissions, meteorology, and atmospheric chemistry (Elliott et al., 2009; Garten, 1992; Hall et al., 2016; Occhipinti, 2008; Russell et al., 1998).

Figure 1 Box (lower quartile, median, upper quartile) and whisker (lower extreme, upper extreme) plot of the distribution of δ¹⁵N values for various NO_x emission sources.

Here we have simulated the emission of ¹⁵NO_xand its mixing in the atmosphere and compared the predicted δ¹⁵N (NO_x, NO${}_{\mathrm{3}}^{-}$) values to observations. The δ¹⁵NO_x values are impacted by three main factors. The first is the inherent variability of the δ¹⁵NO_x emissions in time and space. Secondly, atmospheric processes mix the emitted NO_x, dispersing multiple emission sources within a mixing lifetime relative to the NO_x chemical lifetime (2–7 h), which depends on its concentration and photooxidation chemistry, that also vary in time and by location (Laughner and Cohen, 2019). And thirdly, isotope effects occurring during tropospheric photochemistry may alter the δ¹⁵NO_x emissions as they are transformed from NO_x into NO_y. In this paper, we consider the effects from the first and second considerations, the temporal and spatial variation in NO_x emission, and the impacts from atmospheric transport and deposition processes (source and mixing hypothesis). We accomplish this by incorporating an input dataset of ¹⁵N emissions used in simulations by the CMAQ (Community Multiscale Air Quality) modeling system. In a previous paper we have discussed the impacts of tropospheric photochemistry by incorporating a ¹⁵N chemical mechanism (Fang et al., 2021) into CMAQ. The ultimate goal is to evaluate the accuracy of the NO_x emission inventory using ¹⁵N.

2 Methodology

2.1 Incorporating ¹⁵N into NO_x emission datasets

The EPA trace pollutant emission model SMOKE (Sparse Matrix Operator Kernel Emissions) was used to simulate ¹⁴NO_x and ¹⁵NO_x emissions. ¹⁴NO_x emissions were estimated using the SMOKE model based on the 2002 NEI (National Emission Inventory; United States Environmental Protection Agency, 2014), and ¹⁵N emissions were determined using these ¹⁴NO_x emissions and the corresponding δ¹⁵N values of NO_x sources from previous research (Table 1). Using the definition of δ¹⁵N (‰), ¹⁵NO_x emitted by each SMOKE processing category (area, biogenic, mobile, and point) was calculated by

$$\begin{array}{}\text{(2)}& {}^{\mathrm{15}}{\mathrm{NO}}_x\left(\mathrm{i}\right)={}^{\mathrm{14}}{\mathrm{NO}}_x\left(\mathrm{i}\right)\times {}^{\mathrm{15}}{R}_{{\mathrm{NO}}_x}\left(i\right),\end{array}$$

where ¹⁴NO_x(i) denotes the NO_x emissions for each category (i) obtained from NEI and SMOKE, and ¹⁵R_NOxi is a ¹⁵N emission factor (¹⁵NO_xi $/$ ¹⁴NO_xi) calculated by rearranging Eq. (1):

$$\begin{array}{}\text{(3)}& {}^{\mathrm{15}}{R}_{{\mathrm{NO}}_x}\left(i\right)=\left({\displaystyle \frac{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{NO}}_{x\left(i\right)}}{\mathrm{1000}}}+\mathrm{1}\right)\times \mathrm{0.0036},\end{array}$$

where δ¹⁵NO_x(i) is the δ¹⁵N value of some NO_x source (i= area, biogenic, mobile, and point).

Annual NO_x emissions for 2002 were obtained from the NEI at the county level and were converted into hourly emissions on a 12 km × 12 km grid as previously published (Spak et al., 2007). The modeling domain includes latitudes between 37 and 45^∘ N and longitudes between 98^∘ W and 78^∘ W, which fully covers the Midwestern US (Fig. 2, in yellow). SMOKE categorizes NO_x emissions into four “processing categories”: biogenic, mobile, point, and area (Table 1). The choice of the 2002 version of NEI is, in part, arbitrary. However, to compare the model-predicted δ¹⁵N values with observations, it requires the emission inventory to be relevant to the same timeframe as the δ¹⁵N measurements of the NO_y. The datasets we compare to the model (discussed below) span from 2002 to 2009; thus the 2002 inventory is more relevant than later inventories (2014 onward). The county-level annual ¹⁴NO_x emission for the Midwestern US from NEI was converted to the dataset with hourly ¹⁴NO_x emissions.

Figure 2 The full geographic domain (yellow) and extracted domain (light grayish purple) for the study.

Table 1 The δ¹⁵N values (in ‰) for NO_x emission sources based on SMOKE processing category and NEI sector.

2.1.1 Biogenic ¹⁵NO_x emissions

The NO_x emission from the soil (biogenic) was modeled in SMOKE using standard techniques (details in the Supplement), and the δ¹⁵N values of biogenic NO_x were taken from previous studies. Li and Wang (2008) measured the NO fluxes using dynamic flow chambers for 2 to 13 d after cropland soil was fertilized by either urea (n=9) or ammonium bicarbonate (n=9), and the δ¹⁵NO_x ranged from −48.9 ‰ to −19.8 ‰. Felix and Elliott (2014) used passive samplers to collect NO₂ in a cornfield for 20 d, with low (−30.8 ‰) and high (−26.5 ‰) fertilizer application. Using active samplers, Miller et al. (2018) collected NO₂ between May and June, finding δ¹⁵N ranging from −44.2 ‰ to −14.0 ‰ (n=37); Yu and Elliott (2017) measured −59.8 ‰ to −23.4 ‰ in 15 samples from soil plots in a fallow field 2 weeks after the precipitation. Based on these studies, we adopted an average δ¹⁵N value for NO_x emissions from the soil of −34.3 ‰ (Li and Wang, 2008; Felix and Elliott, 2014; Yu and Elliott, 2017; Miller et al., 2018).

2.1.2 Mobile ¹⁵NO_x emissions

The SMOKE NO_x emission from on-road vehicles used standard methods (details in the Supplement) and used δ¹⁵N values from prior studies. We have excluded studies that infer δ¹⁵ NO_x by measuring plant proxies or passive sampling in the environment (Ammann et al., 1999; Pearson et al., 2000; Savard et al., 2009; Redling et al., 2013; Felix and Elliott, 2014). This is because equilibrium and kinetic isotope effects occur as NO_x reacts in the atmosphere to form NO_y, prior to NO_x deposition. In addition, the role vegetation plays in NO_x removal and atmospheric processes that mix emitted NO_x with the surroundings can also alter the δ¹⁵NO_x. Instead, we estimated the δ¹⁵NO_x emissions from vehicles only using studies that directly measured tailpipe NO_x emissions. Moore (1977) and Heaton (1990) collected tailpipe NO_x spanning −13 ‰ to 2 ‰, with an average of −7.5 ± 4.7 ‰. Neither Heaton nor Moore noted whether these six vehicles were equipped with any catalytic NO_x reduction technology, but it is unlikely since late 1970 and 1980s vehicles were seldom equipped with catalytic NO_x reduction technology. Fibiger (2014) measured five samples of NO_x from diesel engines without SCR (selective catalytic reduction) emitted into a smog chamber; the δ¹⁵N values range from −19.2 ‰ to −16.7 ‰ (±0.97 ‰). The most comprehensive studies on vehicle NO_xδ¹⁵N values are by Walters et al. (2015a, b), who measured gas and diesel vehicles separately, including those with and without three-way catalytic converter (TCC) and SCR technology. They also measured on-road and off-road vehicles separately. The measurements showed that the δ¹⁵NO_x emitted by on-road diesel vehicles ranged from −5 ‰ to 0 ‰, so the average −2.5 ‰ was adopted. The δ¹⁵NO_x values emitted by on-road gasoline vehicles are a function of vehicle travel times, ranging from −6.3 ‰ to 1.8 ‰, with an average of −2.7 ± 0.8 ‰ for the Midwest region. This value is close to the measurements (−8 ‰ to −1 ‰, average −4.7 ± 1.7 ‰) of Miller et al. (2017), who collected NO_x along highways in Pennsylvania and Ohio.

The emission rate of ¹⁵NO_x from the mobile source was determined by Eq. (4) grid by grid, according to the contributions from on-road gasoline vehicles and on-road diesel vehicles, as well as their corresponding δ¹⁵N values. NO_x emissions from off-road vehicles are regarded as area sources in SMOKE, which were processed over each county. In contrast, NO_x emissions from on-road vehicles are regarded as the mobile source in SMOKE, which will be processed along each highway. The δ¹⁵N of on-road gasoline vehicles was based on the average of the vehicle travel time (t) within each region with the same zip code (Walters et al., 2015b).

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{}^{\mathrm{15}}{\mathrm{NO}}_x\left(\text{mobile}\right)& =\left({\displaystyle \frac{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{NO}}_{x\left(\text{on-road gas}\right)}}{\mathrm{1000}}}+\mathrm{1}\right)\\ & \times \mathrm{0.0036}{\times }^{\mathrm{14}}{\mathrm{NO}}_x\left(\text{on-road gas}\right)\\ & +\left({\displaystyle \frac{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{NO}}_{x\left(\text{on-road diesel}\right)}}{\mathrm{1000}}}+\mathrm{1}\right)\\ & \times \mathrm{0.0036}{\times }^{\mathrm{14}}{\mathrm{NO}}_x\left(\text{on-road\hspace{0.17em} diesel}\right),\end{array}\end{array}$$

where ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{NO}}_{x\left(\text{on-road gas}\right)}=-\mathrm{12.35}+\mathrm{3.02}\times \ln (t+\mathrm{0.455})$.

2.1.3 Point source ¹⁵NO_x emissions

NO_x point sources are large anthropogenic NO_x emitters located at a fixed position such as EGUs (electric generating units). Fugitive dust does not significantly contribute to point NO_x emissions, so our inventory focused only on power plants (Houyoux, 2005). Power plants were separated into two different types: EGU and non-EGU (e.g. commercial and industrial combustion). The δ¹⁵N values of NO_x emitted from power plants have been estimated to vary from −19.7 ‰ to 25.6 ‰ (Heaton, 1987, 1990; Snape et al., 2003; Felix et al., 2012; Walters et al., 2015b). We have ignored studies that measured δ¹⁵NO${}_{\mathrm{3}}^{-}$ or δ¹⁵HNO₃ from EGUs (Felix et al., 2015; Savard et al., 2017) and instead, only consider those studies that directly measured δ¹⁵NO_x from stacks. Heaton (1990) collected five samples from the different coal-fired power stations finding NO_x from 6 ‰ to 13 ‰, with a standard deviation of 2.9 ‰. Snape et al. (2003) measured δ¹⁵N values from power plants using three different types of coals values ranging from 2.1 ‰ to 7.2 ‰, with a standard deviation of 1.37 ‰ (n=36). The most comprehensive study on coal-fired power plants NO_x values was by Felix et al. (2012). They measured the δ¹⁵NO_x emission from the coal-fired power stations with and without different emission control technologies. The δ¹⁵NO_x emissions range from 9 ‰ to 25.6 ‰, with an average of 14.2 ± 4.51 ‰ (n=42). The δ¹⁵NO_x values varied when different emission control technologies were used: ranging from 15.5 ‰ to 25.6 ‰, with an average of 19.4 ± 2.28 ‰ (n=16) for SCR (selective catalytic reduction); ranging from 13.6 ‰ to 15.1 ‰, with an average 14.2 ± 0.79 ‰ (n=3) for SNCR (selective noncatalytic reduction); ranging from 9.0 ‰ to 12.6 ‰, with an average 10.7 ± 1.11 ‰ (n=15) for OFA (over-fire air)/LNB (low NO_x burner); and ranging from 9.6 ‰ to 11.7 ‰, with an average 10.5 ± 0.79 ‰ (n=8) for no emission control technology. According to Xing et al. (2013), about half of the coal-fired power plants in the United States are equipped with SCR. Thus, we assume 15 ‰ for the NO_x emissions from coal-fired power plants, which is the average between SCR and other emission control technologies.

The most comprehensive study on natural-gas-fired δ¹⁵NO_x values (Walters et al., 2015b) collected NO_x from a residential natural-gas low-NO_x furnace and the stack of a natural-gas EGU. The measurement showed that the δ¹⁵N values of NO_x emitted by natural-gas power plants averaged −16.5 ± 1.7 ‰, which we used for the NO_x emission from natural-gas power plants. The latitude, longitude, and point source characteristics (EGU and non-EGU, coal-fired or natural-gas-fired, implementation of emission control technology) of each power plant were obtained from the US Energy Information Administration (2017). The power plants were assigned grids by their latitudes and longitudes, and the δ¹⁵N values were assigned to these grids based on their emission characteristics, before determining the emission rate of ¹⁵NO_x from point sources using Eqs. (2) and (3).

2.1.4 Area source ¹⁵NO_x emissions

Area NO_x (details in the Supplement) δ¹⁵N values were based on the assumption that livestock waste and off-road vehicles (utility vehicles for agricultural and residential purposes) accounted for total area sources. Livestock waste δ¹⁵ NO_x values were taken from Felix and Elliott (2014) since it is currently the only study on livestock waste emissions. They placed a passive sampler with ventilation fans in an open-air and closed room in barns of cows and turkeys, respectively. The δ¹⁵NO_x emissions from these measurements range from −29 ‰ to −8.5 ‰. Among these samples, the δ¹⁵NO_x emissions from turkey waste averaged −8.5 ‰, and the δ¹⁵NO_x emissions from cow waste averaged −24.7 ‰. We used −18.8 ‰ as the values of δ¹⁵NO_x emissions from livestock waste, which is the weighted average of the turkey waste and cow waste emissions. We used Walters et al. (2015b) to estimate the δ¹⁵NO_x emissions from the off-road vehicles since it is the latest in-depth study that measured the δ¹⁵NO_x specifically from off-road vehicles that ranged from −15.6 ‰ to −6.2 ‰ and averaged −11.5 ± 2.7 ‰. The measurement showed that the δ¹⁵N values of NO_x emitted by diesel off-road vehicles without SCR ranged from −21.1 ‰ to −16.8 ‰, with an average of −19 ‰ ± 2 ‰, and diesel-powered off-road vehicles with SCR ranged from −9 ‰ to 8.5 ‰, with an average of −2 ‰ ± 8 ‰. We adopted −10.5 ‰ for δ¹⁵N values of NO_x emitted by diesel-powered off-road vehicles, which is the median between the measurement of vehicles with and without SCR.

The emission rate of ¹⁵NO_x from area sources was determined by Eq. (5) grid by grid, according to the contributions from waste, off-road gasoline vehicles, and off-road diesel vehicles, as well as their corresponding δ¹⁵N values based on previous research.

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{}^{\mathrm{15}}{\mathrm{NO}}_x\left(\text{area}\right)& =\left({\displaystyle \frac{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{NO}}_{x\left(\text{waste}\right)}}{\mathrm{1000}}}+\mathrm{1}\right)\\ & \times \mathrm{0.0036}{\times }^{\mathrm{14}}{\mathrm{NO}}_x\left(\text{waste}\right)\\ & +\left({\displaystyle \frac{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{NO}}_{x\left(\text{off-road gas}\right)}}{\mathrm{1000}}}+\mathrm{1}\right)\\ & \times \mathrm{0.0036}{\times }^{\mathrm{14}}{\mathrm{NO}}_x\left(\text{off-road gas}\right)\\ & +\left({\displaystyle \frac{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{NO}}_{x\left(\text{off-road diesel}\right)}}{\mathrm{1000}}}+\mathrm{1}\right)\\ & \times \mathrm{0.0036}{\times }^{\mathrm{14}}{\mathrm{NO}}_x\left(\text{off-road diesel}\right)\end{array}\end{array}$$

The ¹⁵NO_x emission data files of each SMOKE processing category were incorporated into the final dataset based on the δ¹⁵N values from previous research (Table 1) and Eqs. (2)–(5).

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& {\mathit{\delta }}^{\mathrm{15}}{\mathrm{NO}}_{x\left(\text{total}\right)}=\\ & \left({\displaystyle \frac{\frac{{}^{\mathrm{15}}{\mathrm{NO}}_x\left(\text{area}\right){+}^{\mathrm{15}}{\mathrm{NO}}_x\left(\text{biog}\right){+}^{\mathrm{15}}{\mathrm{NO}}_x\left(\text{mobile}\right){+}^{\mathrm{15}}{\mathrm{NO}}_x\left(\text{point}\right)}{{}^{\mathrm{14}}{\mathrm{NO}}_x\left(\text{area}\right){+}^{\mathrm{14}}{\mathrm{NO}}_x\left(\text{biog}\right){+}^{\mathrm{14}}{\mathrm{NO}}_x\left(\text{mobile}\right){+}^{\mathrm{14}}{\mathrm{NO}}_x\left(\text{point}\right)}}{\mathrm{0.0036}}}-\mathrm{1}\right)\\ & \times \mathrm{1000}\end{array}\end{array}$$

2.2 Simulating atmospheric δ¹⁵NO_x in CMAQ

In order to investigate the role of mixing in the spatiotemporal distribution of δ¹⁵NO_x values, CMAQ was used to simulate the meteorological transport effects (advection, eddy diffusion, etc.). In this “emission + transport” scenario, grid-specific δ¹⁵NO_x values emitted are dispersed as NO_x mixes across the regional scale. This dispersion will depend on grid emission strength and mixing vigor and is effectively treating NO_x as a conservative tracer. The simulations used the 2002 National Emission Inventory (NEI), as well as 2002 and 2016 meteorological conditions respectively, to explore how meteorological conditions will impact the atmospheric δ¹⁵NO_x. Simulations covering the full domain and extracted domain were conducted to explore and eliminate potential bias near the domain boundary.

2.2.1 Meteorology input dataset and boundary conditions

To explore the impact of atmospheric processes, the meteorology input datasets for the years 2002 and 2016 were prepared and compared. The CMAQ chemistry-transport model (CCTM) used the NARR (North American Regional Reanalysis) and NAM (North American Mesoscale Forecast System) to convert the weather observations (every 3 h for NARR, every 6 h for NAM analyses) into gridded meteorological elements, such as temperature, wind field, and precipitation, with the horizontal resolution of 12 km and 34 vertical layers, with the thickness increasing with height from 50 m near the surface to 600 m near the 50 mb pressure level. These were used to generate the gridded meteorology files on an hourly basis, using the Weather Research and Forecasting (WRF) model. To maintain consistency between the NO_x emission dataset and the meteorology, the same coordinate system, spatial domain, and grid size used in the SMOKE model were used in the WRF simulation. The WRF outputs were used to prepare the CMAQ-ready meteorology input dataset using CMAQ's MCIP (the Meteorology-Chemistry Interface Processor; see the Supplement for details). In these emission-only simulations, the deposition of NO_x was effectively set to zero. This was accomplished by defining YO =¹⁴NO and YO₂=¹⁴NO₂ (in addition to ZO =¹⁵NO and ZO₂=¹⁵NO₂) and setting their deposition velocities to 0.001 (setting them to zero collapses the simulation). The meteorological fields generated by MCIP were used as inputs for the Initial Conditions Processor (ICON) and Boundary Conditions Processor (BCON) to run CCTM in CMAQ. The ICON program prepares the initial chemical/isotopic concentrations in each of the 3D grid cells for use in the initial time step of the CCTM simulation. The BCON program prepares the chemical/isotopic boundary condition throughout the CCTM simulation. The CMAQ default ICON and BCON for a clean atmosphere were used, which had NO_x < 0.25 ppb. The ¹⁵NO_x values were added to the outputs of ICON and BCON, with the concentration equal to 0.0036[¹⁴NO_x], which assumes δ¹⁵N = 0 at the initial time step and outside the domain of the simulation.

2.2.2 The role of deposition and chemical transformation of NO_x

CMAQ simulated how NO_x removal by photochemical oxidation and deposition alters δ¹⁵NO_x during mixing, transport, and dispersion. This “apparent” conversion of NO_x into NO_y was implemented by enhancing NO_x dry deposition by first magnifying it to 20 times the normal level (14 kg ha⁻¹ yr⁻¹) and testing for the change in NO_x concentration relative to the normal deposition rate. Multiple tuning trials were conducted until the e-folding time (lifetime) of NO_x in the atmosphere across the domain averaged about 1 d. This is a typical average photochemical NO_x lifetime for a combination of urban, suburban, and rural environments (Laughner and Cohen, 2019). This approach is limited since NO_x lifetime varies depending on oxidation capacity, with urban NO_x lifetimes (∼ 2–11 h) being significantly shorter than in rural conditions (Fang et al., 2021). In these simulations, the molecular mass of Y and Z was set equal (14) to ensure no isotope effect was induced by dry deposition, since the equations for dry deposition have a mass term in the diffusion coefficient calculation. These “emission + transport + enhanced NO_x loss” simulations are an attempt to show how “lifetime chemistry” alters δ¹⁵NO_x values by removing NO_x before it can be transported along significant distances. For example, in an emission + transport scenario, NO_x from a high emission power plant could travel across the domain, altering regional δ¹⁵NO_x as it mixes with other grids. In contrast, in the emission + transport + enhanced NO_x loss scenario, most of that NO_x would be removed near the power plant, effectively constricting its δ¹⁵N influence. This has an added advantage in that the deposited δ¹⁵NO_x should be similar to the δ¹⁵NO${}_{\mathrm{3}}^{-}$, which is not being generated in this model. We emphasize that in this model the isotope effects associated with the photochemical transformation of NO_x into HNO₃ (and other higher N oxides) and deposition are ignored and will be addressed in a forthcoming paper.

2.2.3 The simulation over the extracted domain

As mentioned in Sect. 2.2.1, atmospheric δ¹⁵NO_x=0 ‰ for initial conditions and boundary conditions. As a result, a bias may occur along the boundary of the research area and mainly occurs under the following two circumstances – firstly, when the air mass is transported out of the research area (Fig. S1 in the Supplement). Due to the lack of the emission dataset, Canada is considered an “emission-free zone” for this research. As a result, the atmospheric NO_x is diluted, which impacts its δ¹⁵N values, especially for that with extreme δ¹⁵N values (δ¹⁵N $<-$15 ‰ or δ¹⁵N > 5‰). Secondly, the air mass with δ¹⁵NO_x=0 is transported from the emission-free zone into the research area (Fig. S2), and the atmospheric δ¹⁵NO_x is flattened. Therefore, to avoid the bias near the border, the extracted domain that only covers Indiana, Illinois, Ohio, and Kentucky was determined (Fig. 2, in light purple), where the measurements of δ¹⁵N values at NADP (National Atmospheric Deposition Program) sites are available (Mase, 2010; Riha, 2013). The boundary condition for the simulation over the extracted domain is based on the CCTM output of the full-domain simulation (BCON code available on http://www.zenodo.org (last access: 8 December 2020) (https://doi.org/10.5281/zenodo.4311986, Fang, 2020b).

3 Results and discussion

3.1 Simulated spatial variability of NO_x emission rates

We first examine the spatial heterogeneity of the NO_x emission rate for a single time period to illustrate the overall pattern of NO_x emission over the domain (Fig. 3). This is because the δ¹⁵NO_x emission is determined by the fraction of each NO_x source (Eq. 6), which in turn is a function of their emission rates. Since our NO_x emissions are gridded by SMOKE using the NEI, they are, by definition, correct with respect to the NEI. However, a brief discussion of the salient geographic distribution of NO_x emissions and comparisons with other studies is warranted for completeness and as a backdrop for the discussion of NO_x fractions and resulting δ¹⁵N values. We have arbitrarily chosen to sum the NO_x emissions during the April to June time period for this discussion (Fig. 3).

Figure 3 Total NO_x emission in the Midwest between April and June in tonnes of nitrogen per day (t N d⁻¹). High NO_x emissions are associated with major urban areas such as Chicago, Detroit, Minneapolis-St Paul, Kansas City, St. Louis, Indianapolis, and Louisville.

The April to June NO_x emissions ranged from less than 0.01 t N d⁻¹ to more than 15 t N d⁻¹, with the seasonal grid average of 0.904 t N d⁻¹. This average agrees well with estimates in previous studies for the United States, which were between 0.81 and 1.02 t N d⁻¹ (Dignon and Hameed, 1989; Farrell et al., 1999; Selden et al., 1999; Xing et al., 2012). Within 75 % of the geographic domain, the NO_x emissions are relatively low, ranging from between 0 and 0.5 t N d⁻¹ (Fig. S3). Geographically, these grids are in rural areas some distance away from metropolitan areas and highways (Fig. 3). NO_x emissions within about 20 % of the grids are relatively moderate, ranging between 0.5 and 2.0 t N d⁻¹ (Fig. S3). Geographically, these grids are mainly located along major highways and areas with medium population densities (Fig. 3). Urban centers comprise about 5 % of the grids within the geographic domain, and these have high NO_x emissions rates, ranging between 2.0 and 15.0 t N d⁻¹ (Fig. S3). The metropolitan area's average is 5.03 t N d⁻¹, which is nearly 14 times the average emission rate over the rest of the grids within the geographic domain (0.37 t N d⁻¹) due to the high vehicle density associated with high population. The highest emission rates are located within large cities as well as the edge of the east coast metropolitan area (Fig. 3). Summing the NO_x emissions among the grids that encompass these major midwestern cities yields city-level NO_x emission rates that vary from 61.2 t N d⁻¹ (Louisville, KY) to 634.1 t N d⁻¹ (Chicago, IL). These city-level NO_x emission rates (Table S4) agree well with estimates derived from the Ozone Monitoring Instrument (Lu et al., 2015). Grids containing power plants are significant NO_x hotspots within the geographic domain. These account for less than 1 % of the grids, but the NO_x emissions from a single grid that contains a power plant can be as high as 93.4 t N d⁻¹. Geographically, the power plants are mainly located along the Ohio River valley, near other water bodies, and often close to metropolitan areas (Fig. 3). The NO_x emission rates of the major power plants within the Midwest simulated by SMOKE (Table S5) match well with the measurement from the Continuous Emission Monitoring System (CEMS) (de Foy et al., 2015; Duncan et al., 2013; Kim et al., 2009). The geographic distribution of grid-level annual NO_x emission density in our simulation also agrees with the county-level annual NO_x emission density discussed in the 2002 NEI booklet (Fig. S4; United States Environmental Protection Agency, 2018b).

We next examine the spatial heterogeneity of the NO_x source fractions (Fig. 4) for the same time period (April to June). The NO_x fraction (f) is defined as the amount of NO_x from a source category (s) normalized to total NO_x (f_s= NO_x(source) $/$ NO_x(total)). The fraction for anthropogenic NO_x emission is defined as the amount of NO_x from a source category normalized to the sum of NO_x emission from anthropogenic sources (f_s= NO_x(source) $/$ (NO_x(total)-NO_x(biogenic))) Since the δ¹⁵NO_x is determined by the NO_x emission fractions within each grid, it is important to understand where in the domain these fractions differ and why. The area sources, which mainly consist of off-road vehicles, agriculture production, residential combustion, and industrial processes, which are individually too low in magnitude to report as point sources, are fairly uniform in their distribution across the domain.

Figure 4 The geographical distribution of the fraction of NO_x emission from each SMOKE processing category (area, biogenic, mobile, and point) over each grid throughout the Midwest between April and June based on NEI 2002.

The SMOKE simulation shows that the f_s varies significantly across the domain. The average area NO_x emission fraction (f_area) was 0.271 for total NO_x emission and 0.290 for anthropogenic NO_x emission within the Midwest from April to June. The f_area values show a clear spatial variation and range from 0.125 to 0.5 over about 75 % of the grids (Fig. S5). Geographically, the grids with relatively higher f_area are in the rural area away from highways, where agriculture is the most common land use classification. In the states of Wisconsin and Missouri, the f_area is slightly lower due to the higher fraction of NO_x emission from biogenic sources (f_biog). In the states of Pennsylvania and Michigan, the f_area is slightly lower due to the higher fraction of NO_x emission from mobile sources (f_mobile). In addition, the grids with f_area greater than 0.75 are mainly located along the Mississippi River and Ohio River, due to wastewater discharge. The f_bio shows a clear spatial variation and is highest in the western portion of the domain (Fig. 4). The f_bio from April to June is less than 0.5 in more than 90 % of the grids within the geographic domain, with the average of 0.065 (Fig. S5). Geographically, the grids with relatively high f_bio are located in the western regions of the Midwest, away from cities and highway where the density of agricultural acreage and natural vegetation is high. Furthermore, the lowest f_bio values occur in the megacities and along the highways, which agrees well with the land use related to the biogenic emission. The April to June SMOKE simulation shows f_mobile of 0.325 for total NO_x emission and 0.347 for anthropogenic NO_x emission. The f_mobile shows a clear spatial variation, with relatively higher f_mobile located in major metropolitan regions and along the highways, where vehicles have the highest density. The value of f_mobile within the geographic domain distributes evenly on the histogram (Fig. S5). Based on the SMOKE simulation, the fraction of NO_x emission from point sources (f_point) is 0.339 for total NO_x emission and 0.363 for anthropogenic NO_x. The f_point values are obviously highest in grids where the power plants are located, mainly along the Ohio River valley and near other water bodies close to metropolitan areas. The point sources occupy only 4 % of the domain grids, and about one-quarter of the power plants are not on the same grids as highways; thus these grids have a f_point>0.9 NO_x.

3.2 Simulated spatial variability in δ¹⁵NO_x

Using these NO_x emission source fractions, the δ¹⁵NO_x values were simulated, and the spatial heterogeneity of δ¹⁵NO_x for a single time period is discussed. The “emission only” simulation of δ¹⁵NO_x values (at 06:00 UTC on 26 July) ranged from −34.3 ‰ to 14.9 ‰ (Fig. 5a). The majority of the grids have δ¹⁵NO_x values lower than −16.3 ‰, which is due to biogenic NO_x emissions (−34.3‰) in sparsely populated areas where intensive agriculture dominates the land use (Fig. 5a). The δ¹⁵NO_x values for grids containing big cities mainly ranged between −8.75 ‰ and −5 ‰ due to the higher fraction of NO_x emission from on-road vehicles (−2.7‰). Similarly, the δ¹⁵NO_x values for grids ranging between −8.75 ‰ and −5 ‰ resolve major highways. The highest value of δ¹⁵N occurs at the grids, where the coal-fired EGUs (+15‰) and hybrid-fired EGUs are the dominant NO_x source (Fig. 5a).

Figure 5 The δ¹⁵N values of NO_x emission (a emission only scenario) and the δ¹⁵N values of atmospheric NO_x based on NEI 2002 and 2016 meteorology (b emission + transport scenario), at 06:00 UTC on 26 July, are presented by color in each grid. The warmer the color, the higher the δ¹⁵N values of atmospheric NO_x. The feature of the transport inside the white box is shown in Fig. 6.

The effect of atmospheric mixing on the δ¹⁵NO_x spatial distribution was then taken into account by coupling the ¹⁵NO_x emissions to the meteorology simulation. There are significant differences between δ¹⁵NO_x values in the emission only (Fig. 5a) and the emission + transport (Fig. 5b) simulations. While emission only δ¹⁵N pattern shows biogenic NO_x emissions dominating the spatial domain, anthropogenic emissions become dominant over most of the grids in the emission + transport simulations, especially for the grids located around major cities and power plants. In general, as isotopically heavier urban NO_x disperses, the grid average increases from −20.2 ‰ under the emission only scenario to −11.5 ‰ under the emission + transport scenario. Similarly, the NO_x emitted along major highways is transported to the surrounding grids, so that the atmospheric NO_x at the grids around the major highways becomes isotopically heavier relative to the emission only scenario. We define Δδ¹⁵N_transport as the δ¹⁵N difference between emission only and emission + transport scenarios. An example of the Δδ¹⁵N_transport effect can be seen in grids encompassing a plume emanating from southern Illinois' Baldwin Energy Complex (marked with a transparent white box in Fig. 5b) that uses subbituminous coal and bituminous coal as its major energy source. The Δδ¹⁵N_transport in the regions is altered as a function of distance away from the EGU. In this time snapshot (06:00 UTC on 26 July), the northeastwards-propagating plume of NO_x emission from the EGU creates higher δ¹⁵NO_x over 135 km away (Fig. 6).

Figure 6 The Δδ¹⁵N_transport along the plume (colored in dark red to orange inside the white box in Fig. 5b) over the distance from the Baldwin Energy Complex power plant (located at the southwestern border of Illinois).

Figure 7 The geographical distribution of the δ¹⁵N value of total NO_x emissions in each season (winter: January–March; spring: April–June; summer: July–September; fall: October–December) in per mil (‰) throughout the Midwest simulated by SMOKE, based on NEI 2002.

3.3 Seasonal variation in δ¹⁵NO_x

We next examine the temporal heterogeneity of δ¹⁵NO_x values over the domain for emission only and interpret them in terms of changes in NO_x emission fractions as a function of time. The predicted δ¹⁵ NO_x value for total emissions in the Midwest during each season shows a significant temporal variation (Fig. 7). The δ¹⁵NO_x ranged from −35 ‰ to 15 ‰, with the annual average over the Midwest at −6.15 ‰. The maps for different seasons show the obvious changes in δ¹⁵N values over western regions of the Midwest, going from −15 ‰ to −5 ‰ in the spring to −35 ‰ to −15 ‰ in the summer. In order to qualitatively analyze the changes in δ¹⁵NO_x among each season, the values over the grids (Fig. 7) were organized into histograms (Fig. S6). The grids with δ¹⁵NO_x between −35 ‰ and −18 ‰ increase dramatically from less than 10 % during fall (October–December) and winter (January–March) to more than 20 % during spring (April–June) and summer (July–September). The grids with δ¹⁵NO_x between −18 ‰ and −2 ‰ decrease from around 90 % during fall and winter to around 75 % during spring and summer. The significant temporal variation in the δ¹⁵NO_x during different seasons can be quantitatively explained by changing fractions of NO_x emission from the biogenic source in any grid (Fig. S7) using Eq. (6). Unlike other NO_x emission sources, the fraction of NO_x emission from biogenic sources changes significantly among each season within the geographic domain, especially over the rural areas (Fig. S7).

Figure 8 The geographical distribution of the δ¹⁵N value of atmospheric NO_x in each season (winter: January–March; spring: April–June; summer: July–September; fall: October–December) in per mil (‰) throughout the Midwest (with zoomed-in view focusing on Indiana), simulated by CMAQ, based on NEI 2002 and 2016 meteorology.

To qualitatively analyze the changes in the fraction of NO_x emission from biogenic sources among each season, the distributions of the fractions among the same cutoffs as the maps in Fig. S7 were shown in the histograms (Fig. S8). In general, the distribution of the fraction shifts to higher values during spring (April–June) and summer (July–September), indicating the increase of biogenic emissions. During this period, the surface sunlight hours, temperature, and precipitation are relatively higher, and as a result, the canopy coverage of the plants becomes higher, which leads to the increase of the NO_x emission from biogenic sources (Pierce, 2001; Vukovich and Pierce, 2002; Schwede et al., 2005; Pouliot and Pierce, 2009; United States Environmental Protection Agency, 2018a). Besides this, the fertilizer application during this period also increases soil NO_x emissions (Li and Wang, 2008; Felix and Elliott, 2014). As a result, the distribution of δ¹⁵NO_x shifts to lower values during these periods (Fig. 7). The percentage of the grids with the fraction of biogenic emission less than 0.125 decreases dramatically from more than 50 % during fall (October–December) and winter (January–March) to less than 35 % during spring (April–June) and summer (July–September). As the NO_x emission from biogenic source becomes dominant, the percentage of the grids with δ¹⁵NO_x between −35 ‰ and −18 ‰ increases, while the percentage of the grids with values between −18 ‰ and −2 ‰ decreases, which sufficiently explains the trends shown in Fig. 7.

The temporal variation in atmospheric δ¹⁵NO_x is also controlled by the propagation of NO_x emissions, which varies seasonally. The temporal heterogeneity of atmospheric δ¹⁵NO_x under the emission + transport scenario is interpreted in terms of changes in the propagation of NO_x emission as a function of time. The predicted seasonal average δ¹⁵NO_x in the Midwest shows significant variations (Fig. 8). On an annual basis, the emission + transport average δ¹⁵NO_x value was −6.10 ‰, which is similar to the emission only average range, but the range (−19.2 ‰ to 11.6‰) was narrower due to NO_x transport and mixing. The maps for different seasons show the obvious changes in δ¹⁵N values over western regions of the Midwest, from −8.75 ‰ to −5 ‰ in fall and winter to −16.25 ‰ to −12.5 ‰ in spring and summer. The spatial heterogeneity of the δ¹⁵NO_x under the emission + transport scenario (Fig. 8) was compared to that under the emission only scenario (Fig. 7). The difference was defined as Δδ¹⁵N_transport (Fig. S9) and had values that ranged from −21.9 ‰ to 31.2 ‰, with an average of 4.9 ‰. The grids with Δδ¹⁵N_transport between −5 ‰ and 0 ‰ are the urban areas and decrease slightly from about 11 % during fall (October–December) and winter (January–March) to 10 % during spring (April–June) and summer (July–September). The grids with Δδ¹⁵N_transport between 0 ‰ and 5 ‰ are typically in the rural areas that are impacted by the urban NO_x emissions and decrease dramatically from more than 50 % during fall and winter to less than 40 % during spring and summer. The grids with Δδ¹⁵N_transport greater than 5 ‰ are in the rural areas obviously impacted by the urban NO_x emission and increase dramatically from less than 40 % during fall and winter to more than 50 % during spring and summer. Therefore, the impacts from transport and mixing are more obvious during spring and summer (Fig. S10).

The planetary boundary layer (PBL) height is an effective indicator showing whether the pollutants are under synoptic conditions which are favorable for the dispersion, mixing, and transport after being emitted into the atmosphere (Oke, 2002; Shu et al., 2017; Liao et al., 2018; Miao et al., 2019). Comparing the distributions of Δδ¹⁵N_transport values (Fig. S9) with the corresponding PBL height (Fig. S11) for each season, the effects of PBL height on the propagation of the air mass are clearly shown. NO_x emitted by power plants is much higher than the emission rates at the surrounding grids and is a hotspot that impacts the δ¹⁵N values on the surrounding grids. As PBL increases, the emitted NO_x from power plant mixes more effectively with the surrounding grid; thus there are higher δ¹⁵NO_x values along the power plant plume transect. The PBL height changes significantly among each season within the geographic domain, especially over Minnesota, Wisconsin, and Iowa (Fig. S11). The PBL height over these areas increases from less than 250 m above the ground level to more than 625 m a.g.l., during spring and summer, which creates a more favorable synoptic condition for the dispersion, mixing, and transport of the pollutants after being emitted into the atmosphere. As a result, the difference in δ¹⁵N values shifts to higher values, showing the stronger effect of atmospheric processes during spring and summer. In order to qualitatively analyze how PBL height affects the δ¹⁵NO_x along power plant plumes, the domain average PBL height for each month was plotted against δ¹⁵NO_x (Fig. 9a). The δ¹⁵N values along the power plants' plumes and PBL heights over the domain have the same seasonal trend. Interestingly, the “turning point” of the δ¹⁵N values is about 1 month later than the turning point of the PBL heights. The scatter plot (Fig. 9b) shows a strong positive correlation (R²=0.85) between the domain average PBL height and average δ¹⁵N value along the power plants' plumes. The positive correlation between PBL height and propagation of air mass, indicated by the evolution of atmospheric δ¹⁵NO_x in this study, agrees well with the corresponding measurement in megacities in China from previous studies (Shu et al., 2017; Liu et al., 2018; Liao et al., 2018).

Figure 9 The time series plot (a) and the scatter plot (b) of the domain average PBL height (m) and the average δ¹⁵N (‰) value of atmospheric NO_x along the plumes of power plants during each month throughout the Midwest simulated by CMAQ, based on NEI 2002 and 2016 meteorology.

Figure 10 The difference between the δ¹⁵N (‰) value of atmospheric NO_x based on 2016 meteorology and 2002 meteorology (Δδ¹⁵N_2002–2016) during each season (winter: January–March; spring: April–June; summer: July–September; fall: October–December), throughout the Midwest, simulated by CMAQ.

3.4 The simulations based on different meteorology input datasets

The spatial heterogeneity of the δ¹⁵NO_x using 2016 meteorology input dataset was compared to that using 2002 meteorology (Fig. S13). Overall, the simulated δ¹⁵NO_x using 2002 meteorology has the similar geographic distribution and seasonal trend as the 2016 simulation. The difference was defined as Δδ¹⁵N_2002–2016 (Fig. 10) and had values that ranged between −1.25 ‰ and +1.25 ‰ over most of the grids. However, in the western part of the domain, where the biogenic NO_x emission is dominant, the more positive Δδ¹⁵N_2002–2016 values (up to 5 ‰) occur during summer and fall. On the other hand, the more negative Δδ¹⁵N_2002–2016 values (up to −5 ‰) occur along the power plant plume during the same period. The spatial heterogeneity of Δδ¹⁵N_2002–2016 indicates how climate change alters the δ¹⁵NO_x. If we have enough input datasets to generate and compare the seasonal/monthly δ¹⁵NO_x over the past 20+ years, the impacts of anomalies in each meteorology variables could be explored. For the current dataset, a similar comparison between the δ¹⁵NO_x and the corresponding PBL height was conducted for the simulation based on 2002 meteorology (Fig. S14) to show how PBL height changes the evolution of δ¹⁵NO_x. Under the 2002 meteorology, lower PBL height during the winter caused surface δ¹⁵NO_x values along the power plants' plumes to be lower relative to 2016 meteorology. On the other hand, due to the higher PBL height during spring and summer 2002, the δ¹⁵N values decreased through July before ending with relatively higher δ¹⁵N values in December. The scatter plot for the simulation based on 2002 meteorology (Fig. S14b) also shows a strong positive correlation between the domain average PBL height and average δ¹⁵N value along the power plants' plumes, with R²=0.78. The videos of atmospheric δ¹⁵NO_x on an hourly basis throughout the years 2002 and 2016 are available on http://www.zenodo.org (last access: 8 December 2021) (https://doi.org/10.5281/zenodo.4311986, Fang, 2020b).

3.5 The simulation over the extracted domain

Analysis of whether there was a difference between the extracted-domain simulation (Fig. 2) and full-domain simulation was conducted by defining Δδ¹⁵N_{extracted-full} and assessing the bias due to the motion of the air mass across the domain boundary (Fig. S17). The Δδ¹⁵N_{extracted-full} values ranged between −0.25 ‰ and +0.25 ‰ over most of the grids, showing that the difference between the extracted-domain simulation and full-domain simulation of δ¹⁵N values is usually trivial. However, near the southern border of the extracted domain, Δδ¹⁵N_{extracted-full} values are close to +0.75 ‰ (fall and winter) and close to +1.00 ‰ (spring and summer), suggesting the extracted domain may be required for accurate δ¹⁵NO_x simulations

3.6 The role of enhanced NO_x loss

The emission + transport + enhanced NO_x loss simulations significantly alter the δ¹⁵NO_x relative to the “normal deposition” scenarios. Again, the enhanced NO_x loss cases are removing NO_x at rates similar to those by removal via its conversion into HNO₃. Thus, the NO_x deposited is ∼ δ¹⁵NO${}_{\mathrm{3}}^{-}$ (assuming no photochemical isotope effects), and the δ¹⁵NO_x is that in the residual NO_x. The impact of enhanced NO_x loss on the residual NO_x was assessed using Δδ¹⁵N_hi-no, the difference between the δ¹⁵NO_x values under the “enhanced NO_x loss” and “no deposition” scenarios. The Δδ¹⁵N_hi-no range was ±4 ‰ and was especially obvious downwind of the locations with large emission rates, such as power plants or megacities (Fig. 11). This can be explained in a similar fashion to the no deposition scenarios (Fig. S18a), where the dispersion of the isotopically heavier NO_x emission from big cities, major highways, and power plants elevated the δ¹⁵NO_x values in rural areas, and the dispersion of the isotopically lighter biogenic NO_x emission lowered the δ¹⁵NO_x values in the surrounding grids located in the suburb of major cities (Fig. S18b). When enhanced NO_x loss is used, the transport, mixing, and dispersion of local NO_x emissions are restricted to a smaller geographical extent (Fig. S18b), leading to different δ¹⁵NO_x values relative to no deposition. The temporal heterogeneity of Δδ¹⁵N_hi-no over the domain was examined and the impact of enhancing deposition rates of NO_x on the δ¹⁵N of atmospheric NO_x was explored on a seasonal basis (Fig. 12). The seasonal Δδ¹⁵N_hi-no values range from −3.67 ‰ to 5.34 ‰, with an average of 0.51 ‰. The overall pattern of the Δδ¹⁵N_hi-no values shows that due to deposition, the atmospheric NO_x became isotopically lighter over the majority of the grids since EGU and vehicle NO_x is not being transported as far. Conversely, in grids that contain or surround power plants and big cities, the δ¹⁵NO_x increases because it is not as effectively mixing with low δ¹⁵NO_x from nearby grids. The enhanced NO_x loss simulation was used as a proxy to present the isotope effects associated with the “pseudo photochemical transformation” of NO_x into NO_y. The complete isotope effect of tropospheric photochemistry will be addressed in future work, which incorporates ¹⁵N into the chemical mechanism of CMAQ for the simulation.

Figure 11 The Δδ¹⁵N_hi-no values at 18:00 UTC on 25 July.

The δ¹⁵NO_x values of dry deposition (a proxy for δ¹⁵NO${}_{\mathrm{3}}^{-}$) simulated by CMAQ show similar monthly variations and seasonal trends to SMOKE (Fig. S22). The ranges of δ¹⁵NO_x values within each month were narrower, compared to the simulation from SMOKE, with a minimum during February (−8.7 ‰ to −4.4 ‰) and a maximum during August (−11.8 ‰ to −4.2 ‰). The seasonal trend shows low δ¹⁵NO_x values in deposition during summer, with the median around −7.4 ‰ and slightly higher values during winter (median around −6.0 ‰). Therefore, the CMAQ simulation inherits the monthly variations and seasonal trends from SMOKE, while the atmospheric NO_x becomes isotopically heavier, after taking atmospheric mixing and transport into account. As mentioned above, most of the NADP sites are located away from big cities and power plants. Thus, the atmospheric mixing and transport led to the isotopically heavier atmospheric NO_x.

Figure 12 The difference between the δ¹⁵N (‰) values of atmospheric NO_x under the enhanced NO_x loss scenario and no deposition scenario (Δδ¹⁵N_hi-no) during each season (winter: January–March; spring: April–June; summer: July–September; fall: October–December), throughout the Midwest, simulated by CMAQ, based on NEI 2002 and 2016 meteorology.

3.7 Model–observation comparison of δ¹⁵NO_x

In order to evaluate the SMOKE/CMAQ simulations of atmospheric δ¹⁵NO_x, they were compared to two recent studies of δ¹⁵NO_x. The first comparison was relative to rainwater measurements in West Lafayette, IN, from 9 July to 5 August 2016 (Walters et al., 2018). The measured δ¹⁵NO_x values ranged from −33.8 ‰ to 0.2 ‰, with a median of −11.2 ± 8.02 ‰. Under the emission + transport + enhanced NO_x loss scenario using 2016 meteorology, the simulated δ¹⁵NO_x mean (−7.9 ± 2.19 ‰) was 3.3 ‰ less negative than the observations, and the range (−15.9 ‰ to −3.7 ‰) was about half that in the observations (Fig. 13, top, Table S7). The predicted δ¹⁵NO_x was similar regardless of whether 2016 or 2002 meteorology was used but was closer to the measured values compared to the emission only simulations (Fig. 13, top). It is not surprising that the measurements are more negative than the observations because the model does not account for isotope fractionation during the conversion of NO_x into NO_y. Our previous work has shown that the photochemical isotope effect enriches NO_y and depletes NO_x (Fang et al., 2021; Walters and Michalski, 2015), and thus the lower measured δ¹⁵NO_x relative to model is consistent with this isotopic depletion. Our model simulations were also compared to on-road vehicle plume measurement along Midwest highways from 8 to 18 August 2015 (Miller et al., 2017). The box plot also shows more accurate estimation of δ¹⁵N after considering the atmospheric mixing with the emission from surrounding grids (Fig. 13, bottom). Using the emission only scenario, the simulated δ¹⁵NO_x mean was about 3 ‰ more negative than the observations. The predicted δ¹⁵NO_x under the emission + transport + enhanced NO_x loss scenario for these samples along Midwest highways was closer to the measured values, compared to the emission only simulations, regardless of whether 2016 or 2002 meteorology was used. The modeled values are quite close to the observations, suggesting that the photochemical isotope effect is small for these samples. This is not surprising given they were collected on major highways where NO_x concentrations are high and the timescale between collection and emission is small, and thus only a small fraction of emitted NO_x would have been converted to NO_y, minimizing the photochemical isotope effect.

Figure 13 The δ¹⁵NO_x distributions at Lafayette, IN (top), and along Midwest highways (bottom), simulated by SMOKE (a), CMAQ based on 2016 (b) and 2002 meteorology (c), compared with the measured δ¹⁵NO_x (d) (box: lower quartile, median, upper quartile; whisker: lower extreme, upper extreme)

The 30-fold enhanced NO_x loss (see methods) was used to simulate the δ¹⁵N value of NO${}_{\mathrm{3}}^{-}$ deposition (δ¹⁵NO${}_{\mathrm{3}}^{-}$) that was then compared to observations (Fig. 14). As previously noted, rather than explicitly converting NO_x into NO_y via the chemical mechanism in CMAQ, which would require writing an isotope-enabled chemical scheme with appropriate rate constants, we amplified NO_x deposition as a surrogate. This amplification reduced the NO_x lifetime to about 1 d; thus by calculating the δ¹⁵NO_x in the deposition fraction, as opposed to residual NO_x in the atmosphere, we are approximating the δ¹⁵NO${}_{\mathrm{3}}^{-}$ in deposition. The simulated δ¹⁵NO${}_{\mathrm{3}}^{-}$ was compared to NO${}_{\mathrm{3}}^{-}$ collected at NADP sites within Indiana, Illinois, and Ohio in the year 2002 (Table S4). The NEI 2002 and WRF2002 were used for the SMOKE emission model and CMAQ simulations, respectively. The value of deposition was calculated by δ¹⁵NO${}_{\mathrm{3}}^{-}=\mathrm{\Sigma }{f}_{\mathrm{NO}x\mathrm{h}}{\mathit{\delta }}^{\mathrm{15}}$NO_xh, where f_NOxh is the hourly mole fraction of NO_x isotopologue deposited (f_NOxh= NO_xh $/$ NO_xT), and δ¹⁵NO_xh is the δ¹⁵N value of NO_x in deposition. The total NO_x deposited (NO_xT) used to calculate f_NOxh was the amount deposited 5 d prior to the sampling date since the NADP deposition collection integrates the week.

Figure 14 The emission + transport + enhanced NO_x loss CMAQ-predicted δ¹⁵N value of NO_x deposition using NEI 2002 and 2002 meteorology compared to the measured δ¹⁵N of rain NO${}_{\mathrm{3}}^{-}$ at NADP sites within IN, IL, and OH. The photochemical isotope enrichment factor (‰) correction used for each site is noted in the legend.

The δ¹⁵N values of NO_x deposition simulated by CMAQ under the emission + transport + enhanced NO_x loss scenario at each site were compared with the measurements of δ¹⁵N values of NO₃ from prior studies (Mase, 2010; Riha, 2013). While the scatter plot shows a moderate positive correlation between observed and simulated δ¹⁵NO${}_{\mathrm{3}}^{-}$, the simulated value is consistently lighter than the sample δ¹⁵NO${}_{\mathrm{3}}^{-}$ (Fig. 14, top). The magnitude of this negative bias varies among the NADP sites (Fig. S23) and is attributed to isotope fractionation during the conversion NO_x into NO_y, which enriches NO${}_{\mathrm{3}}^{-}$ (Fang et al., 2021; Walters and Michalski, 2015). Globally, this enrichment has been estimated at 3.9 ± 1.8 ‰ (Song et al., 2021). But this enrichment is a function of NO_x, volatile organic compounds (VOCs), and oxidant loading, as well as temperature and photolysis rate (Fang et al., 2021), and is not expected to be the same at each NADP site. After adjusting the simulated δ¹⁵N by raising the values by the average of the difference between sample δ¹⁵N and simulated δ¹⁵N for each site, the scatter plots of sample δ¹⁵N vs. simulated δ¹⁵N fit well into the 1 : 1 line (Fig. 14, bottom). The complete ¹⁵N-incorporated chemical mechanisms will be explored in a future study.

4 Conclusions

The evolution of δ¹⁵N values along the “journey” of atmospheric NO_x were traced, using our ¹⁵N-incorporated SMOKE and CMAQ. The δ¹⁵NO_x under the emission only scenario was simulated by SMOKE, using the NO_x emissions from NEI emission sectors and the corresponding δ¹⁵N values from previous research. The SMOKE simulation indicates that the NO_x emission from biogenic sources is the key driver for the variation of δ¹⁵N, especially among the Midwestern NADP sites. The uncertainties in the δ¹⁵NO_x emission are less than 5 ‰ over the majority of the grids within the Midwest, which were well below the difference among the assigned δ¹⁵NO_x values for different NO_x emission sources (Fig. S24). The δ¹⁵NO_x under the emission + transport scenario was simulated by CMAQ, using the ¹⁵N-incorporated emission input dataset generated from SMOKE, as well as the meteorology input dataset generated from WRF and MCIP. The CMAQ simulation indicates that the PBL height is the key driver for the mixture of anthropogenic and natural NO_x emission, which deepens the gap between δ¹⁵N of atmospheric NO_x and NO_x emission. The δ¹⁵NO_x under the emission + transport + enhanced NO_x loss scenario was simulated by enhancing NO_x deposition in CMAQ simulation, to show how lifetime chemistry alters δ¹⁵NO_x values before it can be transported along significant distances, assuming no isotope fractionation during chemical conversion or deposition.

The simulations under emission only scenario and emission + transport + enhanced NO_x loss scenario were compared to the measurements in West Lafayette, Indiana. The simulated δ¹⁵N agreed well with the seasonal trend and monthly variation. The simulated δ¹⁵NO_x under the emission only scenario was less negative than the corresponding measurements in West Lafayette, IN, taken from July to August 2016. Thus, if we only consider the effects from NO_x emission sources, it is possible that the emissions from soil, livestock waste, off-road vehicles, and natural-gas power plant in West Lafayette, IN, are underestimated, and the emissions from the on-road vehicles and coal-fired power plants in West Lafayette, IN, are possibly overestimated. The simulated δ¹⁵NO_x under the emission + transport + enhanced NO_x loss scenario was about 3 ‰ closer to the corresponding measurements in West Lafayette, IN, compared to the emission only simulations. The simulations under the emission + transport + enhanced NO_x loss scenario were also compared to the measurements of δ¹⁵NO${}_{\mathrm{3}}^{-}$ from NADP sites within Indiana, Illinois, Ohio, and Kentucky. The sample-by-sample comparison shows a moderate positive correlation between observed and simulated δ¹⁵NO${}_{\mathrm{3}}^{-}$, with the negative bias varying among the NADP sites. This bias is attributed to isotope fractionation during the conversion of NO_x into NO_y, affected by different NO_x, VOCs, and oxidant loading, as well as temperature and photolysis rate, at each NADP site. Therefore, future work will explore how tropospheric photochemistry alters δ¹⁵NO_x by incorporating ¹⁵N into the chemical mechanism of CMAQ and comparing the simulation with the corresponding measurements. With the validation of our nitrogen isotopes incorporating CMAQ, the NO_x emission inventories could be effectively evaluated and improved.