Mechanistic representation of soil nitrogen emissions in the Community Multiscale Air Quality (CMAQ) model v 5.1

Soils are important sources of emissions of nitrogen-containing (N-containing) gases such as nitric oxide (NO), nitrous acid (HONO), nitrous oxide (N₂O), and ammonia (NH₃). However, most contemporary air quality models lack a mechanistic representation of the biogeochemical processes that form these gases. They typically use heavily parameterized equations to simulate emissions of NO independently from NH₃ and do not quantify emissions of HONO or N₂O. This study introduces a mechanistic, process-oriented representation of soil emissions of N species (NO, HONO, N₂O, and NH₃) that we have recently implemented in the Community Multiscale Air Quality (CMAQ) model. The mechanistic scheme accounts for biogeochemical processes for soil N transformations such as mineralization, volatilization, nitrification, and denitrification. The rates of these processes are influenced by soil parameters, meteorology, land use, and mineral N availability. We account for spatial heterogeneity in soil conditions and biome types by using a global dataset for soil carbon (C) and N across terrestrial ecosystems to estimate daily mineral N availability in nonagricultural soils, which was not accounted for in earlier parameterizations for soil NO. Our mechanistic scheme also uses daily year-specific fertilizer use estimates from the Environmental Policy Integrated Climate (EPIC v0509) agricultural model. A soil map with sub-grid biome definitions was used to represent conditions over the continental United States. CMAQ modeling for May and July 2011 shows improvement in model performance in simulated NO₂ columns compared to Ozone Monitoring Instrument (OMI) satellite retrievals for regions where soils are the dominant source of NO emissions. We also assess how the new scheme affects model performance for NO_x (NO+NO₂), fine nitrate (NO₃) particulate matter, and ozone observed by various ground-based monitoring networks. Soil NO emissions in the new mechanistic scheme tend to fall between the magnitudes of the previous parametric schemes and display much more spatial heterogeneity. The new mechanistic scheme also accounts for soil HONO, which had been ignored by parametric schemes.

1 Introduction

Global food production and fertilizer use are projected to double in this half-century in order to meet the demand from growing populations (Frink et al., 1999; Tilman et al., 2001). Increasing nitrogen (N) fertilization to meet food demand has been accompanied by increasing soil N emissions across the globe, including in the United States (Davidson et al., 2012). N fertilizer consumption globally increased from 0.9 to 7.4 g N per m⁻² cropland yr⁻¹ between 1961 and 2013, with the US still among the top five N fertilizer users in the world (Lu and Tian, 2017). US N fertilizer use increased from 0.28 to 9.54 g N m⁻² yr⁻¹ during 1940 to 2015. In the past century, hotspots of N fertilizer use have shifted from the southeastern and eastern US to the Midwest and the Great Plains comprising the Corn Belt region (Cao et al., 2017). Recent studies have pointed to soils as a significant source of NO_x emissions, contributing ∼20 % to the total budget globally and larger fractions over heavily fertilized agricultural regions (Jaeglé et al., 2005; Vinken et al., 2014; Wang et al., 2017).

Despite the significance of NO_x emissions generated by soil microbes, policies both globally and for the continental US (CONUS) have focused largely on limiting mobile and point fossil fuel sources of NO_x (Li et al., 2016). Hence, it is incumbent to strategize for the reduction of non-point soil sources of NO_x emissions, especially in agricultural areas. Recent studies have shown higher soil NO_x, even in nonagricultural areas like forests, to significantly impact summertime ozone in CONUS (Hickman et al., 2010; Travis et al., 2016). Consequently, it is increasingly important to estimate both N-fertilizer-induced and nonagricultural NH₃ and NO_x emissions in air quality models.

Soil NO emissions tend to peak in the summertime, when they can contribute 15 %–40 % of the total tropospheric NO₂ column in the continental CONUS (Williams et al., 1992; Hudman et al., 2012; Rasool et al., 2016). Summer is also the peak season for ozone concentrations (Cooper et al., 2014; Strode et al., 2015) and the time when photochemistry is most sensitive to NO_x (Simon et al., 2014). N oxides (NO_x = NO + NO₂) worsen air quality and threaten human health directly and by contributing to the formation of other pollutants. NO_x drives the formation of tropospheric ozone and contributes to a significant fraction of both inorganic and organic particulate matter (PM) (Seinfeld and Pandis, 2012; Wang et al., 2013). Global emissions of NO_x are responsible for one in eight premature deaths worldwide as reported by the World Health Organization (Neira, 2014). The premature deaths are a result of the link of these pollutants to cardiovascular and chronically obstructive pulmonary (COPD) diseases, asthma, cancer, birth defects, and sudden infant death syndrome. These adverse health impacts have been shown to worsen with the rising rate of reactive N emissions from soil N cycling (Kampa and Castanas, 2008; Townsend et al., 2003). NO_x indirectly impacts Earth's radiative balance by modulating concentrations of OH radicals, the dominant oxidant of certain greenhouse gases such as methane (IPCC, 2013; Steinkamp and Lawrence, 2011). Nitrous acid (HONO) upon photolysis releases OH radicals along with NO, driving tropospheric ozone and secondary aerosol formation (Pusede et al., 2012). Soils and agriculture are the leading emitters of N₂O, a potent greenhouse gas (IPCC, 2013).

Ammonia (NH₃) also contributes to a large fraction of airborne fine particulate matter (PM_2.5) (Kwok et al., 2013). Elevated levels of PM_2.5 are linked to various adverse cardiovascular ailments, such as irregular heartbeat and aggravated asthma, that cause premature death (Pope et al., 2009) and contribute to visibility impairment through haze (Wang et al., 2012). NH₃ gaseous emissions also influence the nucleation of new particles (Holmes, 2007). Air quality models such as the Community Multiscale Air Quality (CMAQ) model and GEOS-Chem represent bidirectional NH₃ exchange between the atmosphere and soil–vegetation, analyzed under varied soil, vegetative, and environmental conditions (Cooter et al., 2012; Bash et al., 2013; Zhu et al., 2015).

NO_x, NH₃, HONO, and N₂O are produced from both microbial and physicochemical processes in soil N cycling, predominantly nitrification and denitrification (Medinets et al., 2015; Parton et al., 2001; Pilegaard, 2013; Su et al., 2011). Nitrification is the oxidation of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ to ${\mathrm{NO}}_{\mathrm{3}}^{-}$ whereby intermediate species such as NO and HONO are emitted along with relatively small amounts of N₂O as byproducts. Denitrification is the reduction of soil ${\mathrm{NO}}_{\mathrm{3}}^{-}$; it produces some NO, but predominantly produces N₂O and N₂ (Firestone and Davidson, 1989; Gödde and Conrad, 2000; Laville et al., 2011; Medinets et al., 2015). The fraction of N emitted as NO and HONO relative to N₂O throughout nitrification and denitrification depends on several factors: soil temperature; water filled pore space (WFPS), which in turn depends on soil texture and soil water content; gas diffusivity; and soil pH. HONO is produced during nitrification only and is a source of NO and OH after undergoing photolysis (Butterbach-Bahl et al., 2013; Conrad, 2002; Ludwig et al., 2001; Oswald et al., 2013; Parton et al., 2001; Venterea and Rolston, 2000).

Whether N₂O or N₂ becomes dominant during denitrification depends on the availability of soil ${\mathrm{NO}}_{\mathrm{3}}^{-}$ relative to available carbon (C), WFPS, soil gas diffusivity, and bulk density (i.e., dry weight of soil divided by its volume, indicating soil compaction and/or aeration by O₂). Denitrification rates are quite low even at high soil N concentrations if available soil C is absent. However, the presence of high NO₃ concentrations with sufficient available C is the inhibiting factor for the conversion of N₂O to N₂, keeping N₂O emissions dominant during denitrification (Weier et al., 1993; Del Grosso et al., 2000). Denitrification N₂O emissions are also found to increase with a decrease in soil pH in the range of 4.0 to 8.0 generally (Liu et al., 2010). Fertilizer application and wet and dry deposition add to the soil NH₄ and NO₃ pools, which undergo transformation to emit soil N as intermediates of nitrification and denitrification (Kesik et al., 2006; Liu et al., 2006; Redding et al., 2016; Schindlbacher et al., 2004).

Soil moisture content is the strongest determinant of nitrification and denitrification rates and the relative proportions of various N gases emitted by each. Increasing soil water content due to wetting events such as irrigation and rainfall can stimulate nitrification and denitrification. Nitrification rates peak 2–3 days after wetting, when excess water has drained away and the rate of downward water movement has decreased. Denitrification rates substantially increase and nitrification rates become much slower in wetter soils. This is also influenced by soil texture; for instance, denitrification is favored in poorly drained clay soils and nitrification is favored in freely draining sandy soils (Barton et al., 1999; Parton et al., 2001).

WFPS is a metric that incorporates the above factors. The relative proportions of NO, HONO, and N₂O emitted vary with WFPS. Dry aerobic conditions (WFPS ∼0 %–55 %) are optimal for nitrification, with soil NO dominating soil N gas emissions at WFPS ∼30 %–55 % (Davidson and Verchot, 2000; Parton et al., 2001). HONO emissions have been observed up to WFPS of 40 % and dominate N gas emissions under very dry and acidic soil conditions (Maljanen et al., 2013; Mamtimin et al., 2016; Oswald et al., 2013; Su et al., 2011). Nitrification influences N₂O production within the range of 30 %–70 % WFPS, whereas denitrification dominates N₂O production in wetter soils. Denitrification N₂O is limited by lower WFPS in spite of sufficient available ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and C (Butterbach-Bahl et al., 2013; Del Grosso et al., 2000; Hu et al., 2015; Medinets et al., 2015; Weier et al., 1993). As a result, NO and HONO emissions tend to decrease with increasing water content, whereas N₂O emissions increase subject to available ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and C (Parton et al., 2001; Oswald et al., 2013).

Extended dry periods also suppress soil NO emissions by limiting substrate diffusion while water-stressed nitrifying bacteria remain dormant, allowing N substrate (${\mathrm{NH}}_{\mathrm{4}}^{+}$ or organic N) to accumulate (Davidson, 1992; Jaeglé et al., 2004; Hudman et al., 2010; Scholes et al., 1997). Rewetting of soil by rain reactivates these microbes, enabling them to metabolize accumulated N substrate (Homyak et al., 2016). The resulting NO pulses can be 10–100 times background emission rates and typically last for 1–2 days (Yienger and Levy, 1995; Hudman et al., 2012; Leitner et al., 2017).

Higher soil temperature is critical in increasing NO emission during nitrification under dry conditions. However, N₂O generated in denitrification positively correlates with soil temperature only when WFPS and N substrate availability in soil are not the limiting factors (Machefert et al., 2002; Robertson and Groffman, 2007). Recently, a nearly 38 % increase in NO emitted was observed under dry conditions (∼25 %–35 % WFPS) in California agricultural soils when soil temperatures rose from 30–35 to 35–40 ^∘C (Oikawa et al., 2015). Temperature-dependent soil NO_x emissions may strongly contribute to the sensitivity of ozone to rising temperatures (Romer et al., 2018). Also, some soil NO is converted to NO₂ and deposited to the plant canopy, reducing the amount of NO_x entering the atmosphere (Ludwig et al., 2001).

Mechanistic models of soil N emissions already exist and are used in the Earth science and soil biogeochemical modeling community (Del Grosso et al., 2000; Manzoni and Porporato, 2009; Parton et al., 2001). However, photochemical models like CMAQ have been using a mechanistic approach only for NH₃, while using simpler parametric approaches for NO (Bash et al., 2013; Rasool et al., 2016). Other N oxide emissions like HONO and N₂O are absent from the parametric schemes used in CMAQ (Butterbach-Bahl et al., 2013; Heil et al., 2016; Su et al., 2011). Variability in soil physicochemical properties like pH, temperature, and moisture, along with nutrient availability, strongly control the spatial and temporal trends of soil N compounds (Medinets et al., 2015; Pilegaard, 2013).

The U.S. Environmental Protection Agency (EPA) Air Pollutant Emissions Trends Data show that anthropogenic sources of NO_x (excluding fertilizers) fell by 60 % in the US since 1980, heightening the relative importance of soils. Area sources of NO_x like soils, along with less than expected reduction in off-road anthropogenic sources, are believed to have contributed to a slowdown in US NO_x reductions from 2011–2016 (Jiang et al., 2018). Hence, accurate and consistent representation of soil N is needed to address uncertainties in their estimates.

The parameterized schemes currently implemented in CMAQ for CONUS, like Yienger–Levy (YL) and the Berkeley–Dalhousie Soil NO_x Parameterization (BDSNP), consider only NO expressed as a fraction of total soil N available, without differentiating the fraction of soil N that occurs as organic N, NH₄, or NO₃ (Hudman et al., 2012; Rasool et al., 2016; Yienger and Levy, 1995). Moreover, these parametric schemes classify soil NO emissions as constant factors for different nonagricultural biomes or ecosystems compiled from reported literature and field estimates worldwide (Davidson and Kingerlee, 1997; Steinkamp and Lawrence, 2011; Yienger and Levy, 1995). These emission factors account for the baseline biogenic NO_x emissions in addition to sources from deposition (all biomes) and fertilizer (agricultural land cover only) in the latest BDSNP parameterization (Hudman et al., 2012; Rasool et al., 2016). Despite their limitations, parameterized schemes do distinguish which biomes exhibit low NO emissions (wetlands, tundra, and temperate or boreal forests) from those producing high soil NO (grasslands, tropical savanna or woodland, and agricultural fields) (Kottek et al., 2006; Rasool et al., 2016; Steinkamp and Lawrence, 2011).

The EPA recently coupled CMAQ with the U.S. Department of Agriculture (USDA) Environmental Policy Integrated Climate (EPIC) agroecosystem model. This integrated EPIC–CMAQ framework accounts for a process-based approach for NH₃ by modeling its bidirectional exchange (Nemitz et al., 2001; Cooter et al., 2010; Pleim et al., 2013). The coupled model uses EPIC to simulate fertilizer application rate, timing, and composition. Then, CMAQ estimates the spatial and temporal trends of the soil ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) pool by tracking the ammonium mass balance throughout processes like fertilization, volatilization, deposition, and nitrification (Bash et al., 2013). Using the EPIC-derived soil N pool better represents the seasonal dynamics of fertilizer-induced N emissions across CONUS (Cooter et al., 2012). The coupling with EPIC reduces CMAQ's error and bias in simulating total NH₃ + ${\mathrm{NH}}_{\mathrm{4}}^{+}$ wet deposition flux and ammonium-related aerosol concentrations (Bash et al., 2013). The BDSNP parametric scheme implemented in CMAQ also uses the daily soil N pool from EPIC (Rasool et al., 2016).

Our work builds a new mechanistic approach for modeling soil N emissions in CMAQ based on the DayCENT (Daily version of CENTURY model) biogeochemical scheme (Del Grosso et al., 2000; Parton et al., 2001), integrating nitrification and denitrification mechanistic processes that generate NO, HONO, N₂O, and N₂ under different soil conditions and meteorology. We compare the NO and HONO emissions estimates and associated estimates of tropospheric NO₂ column, ozone, and PM_2.5 with those obtained from CMAQ using the YL and BDSNP parametric schemes. For agricultural biomes, our mechanistic scheme uses daily soil N pools from the same EPIC simulations as in Rasool et al. (2016). Unlike BDSNP, which uses a total weighted soil N, the new mechanistic model tracks different forms of soil N as NH₄, NO₃, and organic N for different soil layers and vegetation types so that nitrification and denitrification can be represented. For nonagricultural biomes, our new mechanistic scheme uses a global soil nutrient dataset in an updated C and N mineralization framework. This enables the model to track the conversion of organic soil N to NH₄ and NO₃ pools on a daily scale for nonagricultural soils.

2 Methodology

2.1 Overview of soil N schemes

Key features of the YL and BDSNP parametric soil NO schemes and our new mechanistic scheme for soil NO, HONO, and N₂O are illustrated in Fig. 1 and Table 1.

Figure 1 Flowchart of the (a) Yienger and Levy (1995) (YL), (b) Berkley–Dalhousie Soil NO_x Parameterization (BDSNP), and (c) mechanistic schemes for soil nitrogen (N) emissions as implemented in CMAQ.

Table 1 Comparison of approaches of the parametric and mechanistic soil N emissions models.

The YL scheme, based on Yienger and Levy (1995), parameterizes soil NO emission ($S_{{\mathrm{NO}}_{\mathrm{YL}}}$, in ng N m⁻² s⁻¹) in Eq. (1) as a function of biome-specific emission factors (A_biome) and soil temperature (T_soil).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{S}_{{\mathrm{NO}}_{\mathrm{YL}}}=f_{\frac{\mathrm{w}}{\mathrm{d}}}\left(A_{\mathrm{biome}(\mathrm{w}/\mathrm{d})},T_{\mathrm{soil}}\right)\\ \text{(1)}& {\displaystyle }& {\displaystyle }P\left(\text{precipitation}\right)\text{CRF}\left(\text{LAI},\text{SAI}\right)\end{array}$$

The emissions factor depends on whether the soil is wet (A_biome(w)) or dry (A_biome(d)), with the wet factor used when rainfall exceeds 1 cm in the prior 2 weeks. For dry soils, YL assumes NO emissions exhibit a small and linear response to increasing soil temperatures. For wet soils, soil NO is zero for frozen conditions, increases linearly from 0 to 10 ^∘C, and increases exponentially from 10 to 30 ^∘C, after which it is constant. In agricultural regions, YL assumes wet conditions throughout the growing season (May–September) and assumes 2.5 % of the fertilizer applied N is emitted as NO, in addition to a baseline NO emissions rate based on grasslands. The pulsing term (P(precipitation)) is applied if precipitation follows at least two dry weeks. The canopy reduction factor (CRF) is set as a function of leaf area index (LAI) and stomatal area index (SAI).

The Biogenic Emissions Inventory System (BEIS v3.61 used in current versions of CMAQ v5.0.2 or higher) estimates NO emissions from soils essentially using the same original YL algorithm as in Eq. (1), with slight updates accounting for soil moisture, crop canopy coverage, and fertilizer application. The YL soil NO algorithm in CMAQ distinguishes between agricultural and nonagricultural land use types (Pouliot and Pierce, 2009). Adjustments due to temperature, precipitation (pulsing), fertilizer application, and canopy uptake are limited to the growing season, assumed as 1 April to 31 October, and are restricted to agricultural areas as defined by the Biogenic Emissions Landuse Database (BELD). Unlike the original YL, the implementation of YL in CMAQ (CMAQ-YL) interpolates between wet and dry conditions based on soil moisture in the top layer (1 cm). In this study, we use the Pleim–Xiu Land Surface Model (PX-LSM) in CMAQ to compute soil temperature (T_soil) and soil moisture (θ_soil).

Agricultural soil NO emissions are based on the baseline grassland NO emission (A_grassland) plus an additional factor (fertilizer(t)) that starts at its peak value during the first month of the growing season and declines linearly to zero at the end of the growing season. The growing season is defined as April–October in CMAQ-YL, rather than being allowed to vary by latitude (original YL) or by a satellite-driven analysis of vegetation (original BDSNP). A summary of the modified YL algorithm is presented below for growing season agricultural emissions (Eq. 2).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{S}_{{\mathrm{NO}}_{\mathrm{CMAQ}-\mathrm{YL}}},\text{agricultural growing season}=\\ {\displaystyle }& {\displaystyle }f\left(A_{\mathrm{grassland}}+\text{fertilizer}\left(t\right),T_{\mathrm{soil}},{\mathit{\theta }}_{\mathrm{soil}}\right)\\ \text{(2)}& {\displaystyle }& {\displaystyle }P\left(\text{precipitation}\right)\mathrm{CRF}\left(\mathrm{LAI},\mathrm{SAI}\right)\end{array}$$

For the nongrowing season or nonagricultural areas throughout the year, soil NO emissions are assumed to depend only on temperature and the base emissions for different biomes (A_biome) as provided in BEIS. CMAQ still uses the base emission for both agricultural and nonagricultural land types with adjustments based solely on air temperature (T_{air, in K}) as done in BEIS (Eq. 3). However, for the sake of simplicity we refer to “CMAQ-YL” merely as “YL”.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{S}_{{\mathrm{NO}}_{\mathrm{CMAQ}-\mathrm{YL}}},\text{nonagricultural or nongrowing season}\\ \text{(3)}& {\displaystyle }& {\displaystyle }=\left(A_{\mathrm{biome}}\right)e^{\left(\mathrm{0.04686}\cdot T_{\mathrm{air}}-\mathrm{14.30579}\right)}\end{array}$$

The original implementation of the BDSNP scheme in CMAQ v5.0.2 was described by Rasool et al. (2016). Here, we update that code for CMAQv5.1, but the formulation remains the same. Soil NO emissions, S_NO, are computed in Eq. (4) as the product of biome-specific emission rates (A_biome(N_avail)) and adjustment factors to represent the influence of ambient conditions. The biome-specific emission rates have background soil NO for 24 MODIS biome types from the literature (Stehfest and Bouwman, 2006; Steinkamp and Lawrence, 2011). Fertilizer and deposition emission rates based on an exponential decay after the input of fertilizer and deposition N are added to background soil NO emission rates for respective biomes. BDSNP accounts for total N from fertilizer and deposition obtained from EPIC. EPIC provides the N available from the crop-specific fertilizer soil N pool in different forms as NH₄, NO₃, and organic N. A final weighted total soil N pool is used by weighting the different N forms by the fraction of each crop type in each modeling grid. The soil temperature response f(T_soil) is an exponential function of temperature (in K). Unlike YL that depends solely on rainfall, BDSNP has a Poisson function g(θ) based on soil moisture (θ) that increases smoothly first until a maximum and then decreases when soil becomes water-saturated. BDSNP also differentiates between wet and dry soil conditions and provides a more detailed representation than YL of pulsing following precipitation and of the CRF (described in Sect. 2.5).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{S}_{{\mathrm{NO}}_{\mathrm{BDSNP}}}=A_{\mathrm{biome}}\left(N_{\mathrm{avail}}\right)f\left(T\right)g\left(\mathit{\theta }\right)P\left(l_{\mathrm{dry}}\right)\\ \text{(4)}& {\displaystyle }& {\displaystyle }\mathrm{CRF}\left(\mathrm{LAI},\text{meteorology},\text{biome}\right)\end{array}$$

Our new mechanistic scheme computes soil emissions of NO, HONO, and N₂O by specifically representing both nitrification and denitrification. Equations (5)–(7) provide an overview of the mechanistic formulation. All functions are described in greater detail in Sect. 2.6.4. In the equations, the pulsing factor P(l_dry) follows the formulation of Rasool et al. (2016). The canopy reduction factor $\text{CRF}\left(\text{LAI},\text{meteorology},\text{biome}\right)$ is described in Sect. 2.5. Briefly, we note that nitrification rates (R_N in Eq. 24, $\mathrm{kg}\text{-}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$) depend on the available NH₄ pool, soil temperature (T_soil), soil moisture (θ_soil), gas diffusivity (Dr), and pH adjustment factors. Meanwhile, denitrification rates (R_D in Eq. (25), kg N ha⁻¹ s⁻¹) depend on the available NO₃ pool, relative availability of NO₃ to C, soil temperature, gas diffusivity, and soil moisture adjustment factors.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{S}_{\mathrm{NO}}=\left({\displaystyle \frac{N_{{\mathrm{NO}}_x}-S_{\mathrm{HONO}}}{\begin{array}{l}+\\ D_{\mathrm{NO}}\end{array}}}\right)\mathrm{CRF}\left(\text{LAI, meteorology, biome}\right)\equiv \\ {\displaystyle }& {\displaystyle }\left(\begin{array}{l}f\left({\mathrm{NH}}_{\mathrm{4}},T_{\mathrm{soil}},{\mathit{\theta }}_{\mathrm{soil}},\mathrm{Dr},\text{pH}\right)P\left(l_{\mathrm{dry}}\right)\\ +\\ f\left({\mathrm{NO}}_{\mathrm{3}}:\mathrm{C},T_{\mathrm{soil}},{\mathit{\theta }}_{\mathrm{soil}},\mathrm{Dr}\right)\end{array}\right)\\ \text{(5)}& {\displaystyle }& {\displaystyle }\mathrm{CRF}\left(\text{LAI, meteorology, biome}\right){\displaystyle }& {\displaystyle }{S}_{\mathrm{HONO}}=\left({\mathrm{HONO}}_f\right)\left(N_{{\mathrm{NO}}_x}\right)\left(f_{\mathrm{SWC}}\right)\\ {\displaystyle }& {\displaystyle }\mathrm{CRF}\left(\text{LAI, meteorology, biome}\right)\equiv \left({\mathrm{HONO}}_f\right)\\ {\displaystyle }& {\displaystyle }\left(f\left({\mathrm{NH}}_{\mathrm{4}},T_{\mathrm{soil}},{\mathit{\theta }}_{\mathrm{soil}},\mathrm{Dr},\mathrm{pH}\right)P\left(l_{\mathrm{dry}}\right)\right)\\ \text{(6)}& {\displaystyle }& {\displaystyle }\left(f_{\mathrm{SWC}}\right)\mathrm{CRF}\left(\text{LAI, meteorology, biome}\right){\displaystyle }& {\displaystyle }{S}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}=\left({\displaystyle \frac{N_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}{\begin{array}{l}+\\ D_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}\end{array}}}\right)\equiv \\ \text{(7)}& {\displaystyle }& {\displaystyle }\left(\begin{array}{l}f\left({\mathrm{NH}}_{\mathrm{4}},T_{\mathrm{soil}},{\mathit{\theta }}_{\mathrm{soil}},\mathrm{Dr},\mathrm{pH}\right)\\ +\\ f\left({\mathrm{NO}}_{\mathrm{3}}:\mathrm{C},T_{\mathrm{soil}},{\mathit{\theta }}_{\mathrm{soil}},\mathrm{Dr}\right)\end{array}\right)\end{array}$$

In all our simulations, soil NH₃ emission is calculated based on the bidirectional exchange scheme (Bash et al., 2013) in CMAQ.

2.2 Biome classification over CONUS

CMAQ uses the National Land Cover Database with 40 classifications (NLCD40; https://www.mrlc.gov/, last access: 22 February 2019) to represent land cover, which is used by the YL parametric scheme. However, Steinkamp and Lawrence (2011) provide soil NO emission factors ($A_{\mathrm{biome}}^{\prime }\left(N_{\mathrm{avail}}\right))$ for only 24 MODIS biomes in the BDSNP parametric scheme. Thus, the initial implementation of BDSNP in CMAQ by Rasool et al. (2016) introduced a mapping between the MODIS 24 and NLCD40 biomes to set an emission factor for each NLCD40 biome type (see Appendix Table A2). Factors were then adjusted using Köppen climate zone classifications (Kottek et al., 2006). Whereas the original implementation of BDSNP by Rasool et al. (2016) treated each grid cell based on its most prevalent biome type, our update of BDSNP for CMAQv5.1 and our mechanistic model use sub-grid biome classification, accounting for the fraction of each biome type in each cell.

The latest Biogenic Emissions Landcover Database version 4 (BELD4), generated using the BELD4 tool in the SA Raster Tools system, is used to represent land cover types consistently across both the Fertilizer Emission Scenario Tool for CMAQ (FEST-C v1.2; https://www.cmascenter.org/fest-c/, last access: 22 February 2019) and the Weather Research and Forecast (WRF) meteorological model (Skamarock et al., 2008) and CMAQ framework. BEIS v3.61 within CMAQ integrates BELD4 with other data sources generated at 1 km resolution to provide fractional crop and vegetation cover. US land use categories are based on the 2011 NLCD40 categories. FEST-C provides tree and crop percentage coverage for 194 tree classes and 42 crops (https://www.cmascenter.org/sa-tools/documentation/4.2/Raster_Users_Guide_4_2.pdf, last access: 22 February 2019). For determining fractional crop cover, the 2011 NLCD–MODIS data were used for Canada and the US in the BELD4 data generation tool of FEST-C. Tree species fractional coverage is based on 2011 Forest Inventory and Analysis (FIA) version 5.1. MODIS satellite products are used where detailed data are unavailable outside of the US.

2.3 N fertilizer

The YL scheme set fertilizer-driven soil NO emissions to be proportional to fertilizer application during a prescribed growing season: May–August for the Northern Hemisphere and November–February for the Southern Hemisphere (Yienger and Levy, 1995) or April–October for CMAQ-YL. Our implementations of both the BDSNP parameterization and mechanistic soil N schemes into CMAQ are designed to enable the use of year- and location-specific fertilizer data with daily resolution. We use FEST-C to incorporate EPIC fertilizer application data into our CMAQ runs. EPIC estimates daily fertilizer application based entirely on simulated idealized plant demand, with N stress and limitations in response to local soil and weather conditions, using linkages with WRF via FEST-C. The FEST-C interface also ensures that EPIC simulations are spatially consistent with CMAQ's CONUS domain and resolution through the Spatial Allocator (SA) Raster Tools system (http://www.cmascenter.org/sa-tools/, last access: 22 February 2019).

Because EPIC covers only the US, outside the US BDSNP uses fertilizer data regridded from Hudman et al. (2012), which scaled Potter et al. (2010) data for fertilizer N from 1994–2001 to global fertilizer levels in 2006. Our mechanistic scheme uses a more recently compiled and speciated soil N and C dataset for non-US agricultural regions, regridded from Xu et al. (2015).

2.4 N deposition

N deposition serves as a significant addition to the soil mineral N (inorganic $\mathrm{N}:{\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$) pool and hence influences soil N emissions. The YL scheme does not explicitly represent N deposition but instead sets soil emissions based on biome type. In our implementation of both the updated BDSNP and new mechanistic soil N schemes, hourly wet and dry deposition rates for both reduced and oxidized forms of N, computed within the CMAQ simulation, are added to the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ soil pools.

2.5 Canopy reduction factor (CRF)

CRF is used to calculate above-canopy NO and HONO, assuming that some fraction of each is converted to NO₂ and absorbed by leaves. Earlier global-scale GEOS-Chem simulations with BDSNP had a monthly averaged CRF that reduced total soil NO_x by an average of 16 % (Hudman et al., 2012).

The original YL soil NO scheme (Yienger and Levy, 1995) and the in-line BEIS in CMAQ set CRF as a function of LAI and SAI. Recently, implementations of BDSNP in CMAQ and GEOS-Chem implemented CRF as a function of wind speed, turbulence, and canopy structure (Geddes et al., 2016; Rasool et al., 2016; Wang et al., 1998).

Here, we compute CRF using equations from Wang et al. (1998) for both the BDSNP and new mechanistic scheme using spatially and temporally variable land-surface parameters: surface (2 m) temperature, solar radiation (W m⁻²), surface pressure, snow cover, wind speed (v_wind), cloud fraction, canopy structure, vegetation coverage (LAI and canopy resistances), gas diffusivity, and deposition coefficients. The final reduction factor (CRF(LAI, meteorology, biome)) for primary biogenic soil NO emissions is based on two main factors: bulk stomatal resistance (R_Bulk) and the land-use-specific ventilation velocity of NO (v_vent,NO), calculated based on the parameters mentioned above (Eq. 8).

$$\begin{array}{}\text{(8)}& {\displaystyle }\mathrm{CRF}\left(\text{LAI, meteorology, biome}\right)={\displaystyle \frac{R_{\mathrm{Bulk}}}{R_{\mathrm{Bulk}}+v_{\mathrm{vent},\mathrm{NO}}}}\end{array}$$

The ventilation velocity of NO (v_vent,NO) is calculated by adjusting a normalized day- and night-specific velocity from Wang et al. (1998): 10⁻² and $\mathrm{0.2}\times {\mathrm{10}}^{-\mathrm{2}}$ m s⁻¹, respectively. The adjustments are based on biome-specific LAI and canopy wind extinction coefficients (C_Biome). C_{tropical rainforest} is the canopy wind extinction coefficient for tropical rainforests, the biome on which most canopy uptake studies for NO_x are based (Eq. 9).

$$\begin{array}{ll}{\displaystyle }{v}_{\mathrm{vent},\mathrm{NO}}=& {\displaystyle }{v}_{\mathrm{vent},{\mathrm{NO}}_{\mathrm{day}/\mathrm{night}}}\sqrt{{\left({\displaystyle \frac{v_{\mathrm{wind}}}{\mathrm{3}}}\right)}^{\mathrm{2}}\left({\displaystyle \frac{\mathrm{7}}{\text{LAI}}}\right)}\\ \text{(9)}& {\displaystyle }& {\displaystyle }\left({\displaystyle \frac{C_{\text{tropical rainforest}}}{C_{\mathrm{Biome}}}}\right)\end{array}$$

R_Bulk is a combination of various canopy resistances in series and parallel: internal stomatal resistance, cuticle resistance, and aerodynamic resistance, which have biome-specific normalized values for the MODIS 24 biomes also available in the dry deposition scheme of CMAQ. These normalized values of individual resistances are subsequently adjusted and dependent on multiple conditions for solar radiation, surface temperature, pressure, deposition coefficients, and the molecular diffusivity of NO₂ in air. The calculation of R_Bulk based on Wang et al. (1998) has been documented and shared in the open-source BDSNP code repository (canopy_nox_mod.F) for the purpose of reproducibility (available at https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=_1351, last access: 22 February 2019).

2.6 Detailed description of the mechanistic soil N scheme

2.6.1 Overview

Our new mechanistic soil N model tracks the NH₄, NO₃, and organic C and N pools in soil separately, in contrast to the total N pool of BDSNP, and estimates NO, HONO, and N₂O rather than just NO (Fig. 2). It uses DayCENT to represent both nitrification and denitrification. For agricultural biomes, we use speciated N and C pools from EPIC to drive DayCENT. For nonagricultural biomes, we use a C–N mineralization framework (Manzoni and Porporato, 2009) to estimate the inorganic N and C pools for DayCENT.

Figure 2 Schematic for N transformation to estimate soil pools of ammonium (NH₄) and nitrate (NO₃) and the resultant nitrification and denitrification N emissions in the mechanistic model.

One of the advantages of using DayCENT is its ability to simulate all types of terrestrial ecosystems. DayCENT is one of the only biogeochemical models that not only provides a process-based representation of soil N emissions, but has also been calibrated and validated across an array of conditions for crop productivity, soil C, soil temperature and water content, N₂O, and soil ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (Necpálová et al., 2015). Hence, mechanistic models like DayCENT yield more reliable results by applying validated controls of soil properties like soil temperature and moisture, which are the key process controls to nitrification and denitrification. More recent mechanistic models like DNDC, MicNit, ECOSYS, and COUPMODEL are quite similar to DayCENT in their representation of the nitrification and denitrification process. However, these models have not been as widely evaluated and impose greater computational costs (Butterbach-Bahl et al., 2013). DayCENT also enhances consistency in our mechanistic model by utilizing the same C–N mineralization scheme (taken from the CENTURY model; Parton et al., 2001) that is used in EPIC.

Most stand-alone applications of DayCENT and other mechanistic models have focused on the biogeochemical, climate, and agricultural impacts of soil emissions. Our linkage of DayCENT with CMAQ provides an opportunity for the first time to estimate emissions of multiple soil N species through a process-based approach and then assess their impact on atmospheric chemistry in a regional photochemical model.

2.6.2 Agricultural regions

In agricultural regions, we use EPIC to derive organic N, NH₄, NO₃, and C pools updated on a daily scale. EPIC follows the same approach used in the CENTURY model (Parton et al., 1994), but uses an updated crop growth model and better represents the effects of sorption on soil water content that affect leaching losses and the surface-to-subsurface flow of N. In contrast, CENTURY used monthly water leached below 30 cm of soil depth, annual precipitation, and the silt and clay content of soil (Izaurralde et al., 2006).

In EPIC, organic N residues added to the agricultural soil surface or belowground from plant or crop residues, roots, fertilizer, deposition, and manure are split into two broad compartments: microbial or active biomass and slow or passive humus. Slow or passive humus is essentially recalcitrant and nonliving in nature with very slow turnover rates ranging from centuries to even thousands of years and makes up most of the organic matter. N uptake by soil microbes from organic matter, also called “microbial biomass” or “microbial–active N”, is the living portion of the soil organic matter, excluding plant roots and soil animals larger than $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$ µm³. Although microbial biomass constitutes a small portion of organic matter (∼2 %), it is central in microbial activity: in other words, the conversion of organic N to inorganic N (Cameron et al., 2013; Manzoni and Porporato, 2009). The transformation rate of organic N to microbial N is controlled by the relative C and N content in microbial biomass, soil temperature and water content, soil silt and clay content, organic residue composition enhanced by tillage in agricultural soil, bulk density, oxygen content, and inorganic N availability. Microbial N has quicker turnover times ranging from days to weeks compared to hundreds of years for slow or passive organic matter (Izaurralde et al., 2006; Schimel and Weintraub, 2003). Hence, microbial biomass is the main clearinghouse and driver of C and N cycling in EPIC. Whether net mineralization of organic N to ${\mathrm{NH}}_{\mathrm{4}}^{+}$ occurs or net immobilization of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ to microbial N depends strongly on the relative C and N contents in microbial biomass. Higher N content supports net mineralization, whereas higher C content supports net immobilization. C and N can also be leached or lost in gaseous forms (Izaurralde et al., 2012).

We then estimate gaseous N emissions by using the organic N, NH₄, NO₃, and C pools provided from EPIC/FEST-C along with relevant soil properties for agricultural biomes from the DayCENT nitrification and denitrification sub-model, as described in Sect. 2.6.4 and illustrated in Fig. 2.

2.6.3 Nonagricultural regions

We adapt the framework for linked C and N cycling from Schimel and Weintraub (2003) for nonagricultural regions, where EPIC is not applicable. This framework accounts for the mineralization of organic N by considering which element is limiting based on the relative C-to-N content in microbial biomass. If N is in excess, then the mineralization of organic N producing ${\mathrm{NH}}_{\mathrm{4}}^{+}$ is favored. If C is in excess, it results in overflow metabolism that results in elevated C respiration rates not associated with microbial growth. The resultant inorganic N and C respiration rates are then applied on a temporal and spatial scale consistent with those for the EPIC agricultural pool.

To ensure mass balance, enzyme production (Eqs. 11–13) and recycling mechanisms (Eqs. 14–15) to replenish microbial biomass C are crucial. Similarly, net immobilization is assumed as was done in EPIC when we approach C-saturated conditions with time to replenish microbial N. Without such mechanisms, there is a danger to always incorrectly predict the N- or C-limited state for microbes. Also, some proportion of the microbial biomass is utilized for the maintenance of living cells (only C demand) (Eq. 14), while the rest accounts for decay and regrowth (both C and N demands) (Eqs. 16–17, 18–19) (Schimel and Weintraub, 2003; Manzoni and Porporato, 2009). Fractions of C and N in dying microbial biomass are recycled into the available microbial C and N pools. Schimel and Weintraub (2003) provide values for parameters that quantify these growth and decay processes: fraction of biome C to exoenzymes (K_e) =0.05; microbial maintenance rate (K_m) =0.01 d⁻¹; substrate use efficiency (SUE) =0.5; proportion of microbial biomass that dies per day (K_t) =0.012 d⁻¹; proportion of microbial biomass (C or N) for microbial use (K_r) =0.85.

$$\begin{array}{ll}\text{(10)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{m}}\left(\text{respiration from maintenance}\right)={\mathrm{K}}_{\mathrm{m}}\left(\text{SMC}\right){\displaystyle }& {\displaystyle }{R}_{\mathrm{e}}\left(\text{respiration from enzyme production}\right)=\\ \text{(11)}& {\displaystyle }& {\displaystyle }\left(\right(\mathrm{1}-\mathrm{SUE}\left)\right({\mathrm{EP}}_{\mathrm{C}})/\mathrm{SUE}){\displaystyle }& {\displaystyle }{\mathrm{EP}}_{\mathrm{C}}\left(\text{enzyme production as, C loss/sink}\right)=\\ \text{(12)}& {\displaystyle }& {\displaystyle }{\mathrm{K}}_{\mathrm{e}}\left(\mathrm{SMC}\right){\displaystyle }& {\displaystyle }{\mathrm{EP}}_{\mathrm{N}}\left(\text{enzyme production as N loss/sink}\right)={\mathrm{EP}}_{\mathrm{C}}/\mathrm{3}\\ \text{(13)}& {\displaystyle }& {\displaystyle }\left(\text{where 3 is the approximate}\mathrm{C}:\mathrm{N}\text{ratio for protein}\right){\displaystyle }& {\displaystyle }{\mathrm{CY}}_{\mathrm{C}}\left(\text{recycle from C microbial biomass}\right)=\\ \text{(14)}& {\displaystyle }& {\displaystyle }{\mathrm{K}}_t{\mathrm{K}}_r\left(\mathrm{SMC}\right){\displaystyle }& {\displaystyle }{\mathrm{CY}}_{\mathrm{N}}\left(\text{recycle from N microbial biomass}\right)=\\ \text{(15)}& {\displaystyle }& {\displaystyle }{\mathrm{CY}}_{\mathrm{C}}/C_{\mathrm{m}}:N_{\mathrm{m}}\text{(16)}& {\displaystyle }& {\displaystyle }{H}_{\mathrm{C}}\left(\text{C death/decay}\right)={\mathrm{K}}_t\left(\mathrm{1}-K_r\right)\left(\mathrm{SMC}\right)\text{(17)}& {\displaystyle }& {\displaystyle }{H}_{\mathrm{N}}\left(\text{N death/decay}\right)=H_{\mathrm{C}}/{\mathrm{C}}_{\mathrm{m}}:{\mathrm{N}}_{\mathrm{m}}{\displaystyle }& {\displaystyle }\text{If C is limited or N in excess:}\\ {\displaystyle }& {\displaystyle }\mathrm{SMC}<R_{\mathrm{m}}+\left({\mathrm{EP}}_{\mathrm{C}}/\mathrm{SUE}\right)\\ \text{(18)}& {\displaystyle }& {\displaystyle }+\left(\left(\mathrm{SMN}-{\mathrm{EP}}_{\mathrm{N}}\right)\left({\mathrm{C}}_{\mathrm{m}}:{\mathrm{N}}_{\mathrm{m}}/\mathrm{SUE}\right)\right){\displaystyle }& {\displaystyle }{R}_{\mathrm{g}}\left(\text{respiration from growth, C limited}\right)=\\ \text{(19)}& {\displaystyle }& {\displaystyle }\left(\mathrm{1}-\mathrm{SUE}\right)\left(\mathrm{SMC}-\left({\mathrm{EP}}_{\mathrm{C}}/\mathrm{SUE}\right)-R_{\mathrm{m}}\right)\text{(20)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{O}}\left(\text{respiration from overflow mechanism}\right)=\mathrm{0}{\displaystyle }& {\displaystyle }{\mathrm{NH}}_{\mathrm{4}}\left(\text{from net mineralization after mass balance}\right)=\\ {\displaystyle }& {\displaystyle }(\mathrm{SMN}-{\mathrm{EP}}_{\mathrm{N}}-(\left(\mathrm{SMC}-\left({\mathrm{EP}}_{\mathrm{C}}/\mathrm{SUE}\right)-R_{\mathrm{m}}\right)\\ \text{(21)}& {\displaystyle }& {\displaystyle }\left(\mathrm{SUE}/{\mathrm{C}}_{\mathrm{m}}:{\mathrm{N}}_{\mathrm{m}}\right)\left)\right)\end{array}$$

We represent spatial heterogeneity in soil C and N by using the Schimel and Weintraub (2003) algorithm with sub-grid land use fractions from NLCD40 to estimate the different parameters for specific nonagricultural biomes in Eqs. (10)–(20). That allows us to account for inter-biome variability in soil properties and organic and/or microbial biomass.

Mineralized N pools generated as ${\mathrm{NH}}_{\mathrm{4}}^{+}$ in this framework are calculated eventually as a function of microbial biomass and the aforementioned parameters driving the net mineralization (Eqs. 18 and 21).

We map a global organic C and N pool dataset (Xu et al., 2015) onto our CONUS domain using biome-specific fractions from 12 different biome types for the conversion of these organic pools into microbial biomass pools (Xu et al., 2013). We map these 12 broader biome types to the 24 MODIS biome types with the mapping shown in Table A1. To ensure consistency with the sub-grid biome fractions for the 40 NLCD biome types (Sect. 2.2), we map the MODIS 24 biome-specific microbial ∕ organic C and N fractions to NLCD 40 (Cmic_biome and Nmic_biome; biome represents the 40 NLCD categories) with the mappings shown in Tables A2 and A3. We calculate area-weighted microbial C and N pools (SMC and SMN) using Cmic_biome and Nmic_biome that account for the inter-biome variability in the availability of soil microbial biomass. Also, spatial heterogeneity in terms of vertical stratification is crucial as emission losses from N cycling primarily happen in the top 30 cm layer. Hence, we incorporate the Xu et al. (2015) data for the top 30 cm for the organic nutrient pool and microbial C:N ratio (C_m:N_m) along with other soil properties such as soil pH, θ_soil, and T_soil. This framework (Fig. 2) enables us to estimate soil NH₄, NO₃, and C pools from area-weighted microbial biomass as consistently as possible with the pools that EPIC provides in agricultural regions.

2.6.4 DayCENT representation of soil N emissions

The final part of the mechanistic framework is formed by using a nitrification and denitrification N emissions sub-model adapted from DayCENT along with nitrification and denitrification rate calculations adapted from EPIC. Nitrification and denitrification rates are adapted from EPIC to maintain consistency with the NH₃ bidirectional scheme in CMAQ, which uses the same. It should be noted that the coupled C–N decomposition module in the EPIC terrestrial ecosystem model is similar to that of DayCENT (Izaurralde et al., 2012, 2017; Gaillard et al., 2018). EPIC-simulated agricultural NH₄ and NO₃ soil pools are generated as described in Sect. 2.6.2, whereas the nonagricultural NH₄ and NO₃ soil pools are calculated by using the methods described in Sect. 2.6.3 (Eqs. 22–23). NH₄ and NO₃ soil pools drive nitrification and denitrification as shown in Eqs. (24)–(25). Variability in terms of the soil conditions influencing N emissions in nitrification and denitrification is introduced through the rates at which NH₄ is nitrified (R_N) and NO₃ is denitrified (R_D) (Eqs. 24–25).

The nitrification rate (K_N) (Eq. 26) is estimated based on regulators from the soil water content, soil pH, and soil temperature (T_soil), following the approach of Williams et al. (2008), consistent with the bidirectional NH₃ scheme in CMAQ (Bash et al., 2013). The nitrification soil temperature regulator (f_T) accounts for frozen soil with no evasive N fluxes (Eq. 27). The nitrification soil water content regulator (f_SW) accounts for soil water content at the wilting point and field capacity (Eqs. 28–29). The regulator terms f_T and f_SW both get their dependent variables from land-surface outputs derived from the Meteorology–Chemistry Interface Processor (MCIP) (Otte and Pleim, 2010). However, the nitrification soil pH regulator (f_pH) takes soil pH for agriculture soil from EPIC and for nonagricultural soil from a separate global dataset (Xu et al., 2015), available at both 0.01 and 1 m depths to maintain consistency with MCIP (Eq. 30). The denitrification rate (K_D) (Eq. 31) is regulated by soil temperature (Eq. 34), with WFPS (Eq. 33) acting as a proxy for O₂ availability and soil moisture (θ_soil ) and the relative availability of NO₃ and C (Eq. 32) determining N₂O or N₂ emissions during denitrification (Williams et al., 2008). Note that Eqs. (26) and (31) set upper limits for K_N and K_D, respectively.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{NO}}_{\mathrm{3}}(\mathrm{kg}\text{-}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}},\text{after nitrification})=\\ \text{(22)}& {\displaystyle }& {\displaystyle }{\mathrm{NH}}_{\mathrm{4}}(\mathrm{1.0}-e^{-\left(K_{\mathrm{N}}\mathrm{d}t\right)}){\displaystyle }& {\displaystyle }{\mathrm{NH}}_{\mathrm{4}}(\mathrm{kg}\text{-}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}},\text{after nitrification})=\\ \text{(23)}& {\displaystyle }& {\displaystyle }{\mathrm{NH}}_{\mathrm{4}}{e}^{-\left(K_{\mathrm{N}}\mathrm{d}t\right)}\text{(24)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{N}}\left(\mathrm{kg}\text{-}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}\right)={\mathrm{NH}}_{\mathrm{4}}(\mathrm{1.0}-e^{-\left(K_{\mathrm{N}}\mathrm{d}t\right)})/\mathrm{d}t\text{(25)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{D}}\left(\mathrm{kg}\text{-}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}\right)={\mathrm{NO}}_{\mathrm{3}}(\mathrm{1.0}-e^{-\left(K_{\mathrm{D}}\mathrm{d}t\right)})/\mathrm{d}t\text{(26)}& {\displaystyle }& {\displaystyle }{K}_{\mathrm{N}}\left({\mathrm{s}}^{-\mathrm{1}}\right)=min\left(\mathrm{0.69},\left(f_{\mathrm{T}}\right)\left(f_{\mathrm{SW}}\right)\left(f_{\mathrm{pH}}\right)\right){\displaystyle }& {\displaystyle }{f}_{\mathrm{T}}\left(\text{nitrification soil temperature regulator}\right)=\\ \text{(27)}& {\displaystyle }& {\displaystyle }max\left(\mathrm{0.041}\left(T_{\mathrm{soil}}-\mathrm{278.15}\right),\mathrm{0.0}\right){\displaystyle }& {\displaystyle }{f}_{\mathrm{SW}}\left(\text{nitrification soil water content regulator}\right)=\\ \text{(28)}& {\displaystyle }& {\displaystyle }\left\{\begin{array}{l}\mathrm{0.1},\text{if}\left({\mathit{\theta }}_{\mathrm{soil}}\le \text{wilting point}\right)\\ max\left(\mathrm{0.1},\mathrm{0.1}+\mathrm{0.9}\sqrt{{\displaystyle \frac{\left({\mathit{\theta }}_{\mathrm{soil}}-\text{wilting point}\right)}{\left(\text{field capacity-wilting point}\right)}}},\right.\\ \left.{\displaystyle \frac{\left({\mathit{\theta }}_{\mathrm{soil}}-\text{wilting point}\right)}{\mathrm{0.25}\left(\text{field capacity–wilting point}\right)}}\right),\\ \text{if}\left(\text{wg25}>{\mathit{\theta }}_{\mathrm{soil}}>\text{wilting point}\right)\\ \mathrm{1.0},\text{if}\left(\text{field capacity}>{\mathit{\theta }}_{\mathrm{soil}}\ge \text{wg25}\right)\\ max\left(\mathrm{0.1},\mathrm{1.0}-{\displaystyle \frac{\left({\mathit{\theta }}_{\mathrm{soil}}-\text{field capacity}\right)}{\left({\mathit{\theta }}_{\mathrm{soil}}\left(\text{at saturation}\right)-\text{field capacity}\right)}}\right),\\ \text{if}\left({\mathit{\theta }}_{\mathrm{soil}}>\text{field capacity}\right)\\ \end{array}\right.\text{(29)}& {\displaystyle }& {\displaystyle }\mathrm{\text{wg25}}=\text{wilting point}+\mathrm{0.25}\left(\text{field capacity–wilting point}\right){\displaystyle }& {\displaystyle }{f}_{\mathrm{pH}}\left(\text{nitrification soil pH regulator}\right)=\\ \text{(30)}& {\displaystyle }& {\displaystyle }\left\{\begin{array}{l}\mathrm{0.307}\left(\text{pH}\right)-\mathrm{1.269},\text{acidic soil}\left(\text{pH}<\mathrm{7}\right)\\ \mathrm{1.0},\text{neutral soil}\left(\mathrm{7.4}>\text{pH}\ge \mathrm{7}\right)\\ \mathrm{5.367}-\mathrm{0.599}\left(\text{pH}\right),\text{alkaline soil}\left(\text{pH}\ge \mathrm{7.4}\right),\end{array}\right.\text{(31)}& {\displaystyle }& {\displaystyle }{K}_{\mathrm{D}}\left({\mathrm{s}}^{-\mathrm{1}}\right)=min\left(\mathrm{0.01},f(\text{WFPS},T_{\mathrm{soil}},{\mathrm{NO}}_{\mathrm{3}}:\mathrm{C})\right){\displaystyle }& {\displaystyle }f\left(\text{WFPS},T_{\mathrm{soil}},{\mathrm{NO}}_{\mathrm{3}}:\mathrm{C}\right),\text{denitrification regulators}=\\ \text{(32)}& {\displaystyle }& {\displaystyle }\left(f_{\mathrm{T},\mathrm{D}}\right)\left(f_{\mathrm{WFPS},\mathrm{D}}\right)\left({\displaystyle \frac{\left(\mathrm{1.4}\left(\text{labile C}\right)\left({\mathrm{NO}}_{\mathrm{3}}\right)\right)}{\left(\left(\text{labile C}+\mathrm{17}\right)\left({\mathrm{NO}}_{\mathrm{3}}+\mathrm{83}\right)\right)}}\right)\text{(33)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{WFPS},\mathrm{D}}=min\left(\mathrm{1.0},{\displaystyle \frac{\mathrm{4.82}}{{\mathrm{14}}^{(\mathrm{16}/({\mathrm{12}}^{\mathrm{1.39}\left(\text{WFPS}\right))})}}}\right)\text{(34)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{T},\mathrm{D}}=min\left(\mathrm{1.0},e^{\left(\mathrm{308.56}\left(\frac{\mathrm{1}}{\mathrm{68.02}}-\frac{\mathrm{1}}{T_{\mathrm{soil}}\left(\text{in K}\right)-\mathrm{227.13}}\right)\right)}\right)\end{array}$$

DayCENT partitions N emissions as NO_x and N₂O based on relative gas diffusivity in soil compared to air (Dr) (Eq. 35). Dr is calculated based on the algorithm from Moldrup et al. (2004), which accounts for soil water content, soil air porosity, and soil type. Dr, and hence the ratio of NO_x to N₂O emissions ($r_{{\mathrm{NO}}_x/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$) being a function of Dr, also accounts for soil texture by quantifying pore space, which is highest in coarse soil (Parton et al., 2001; Moldrup et al., 2004). DayCENT assumes 2 % of nitrified N (R_N) is lost as N₂O (Eq. 36). $r_{{\mathrm{NO}}_x/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$ is the ratio of NO_x (both NO and HONO, which photolyze rapidly to NO) to N₂O, in which emissions are expressed on a g N h⁻¹ basis. These emissions are susceptible to pulsing after rewetting of soil in arid or semiarid conditions (P(l_dry)), as explained in Sect. 2.1 (Eq. 37). Denitrification NO is also calculated using the overall $r_{{\mathrm{NO}}_x/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$ ratio (Eq. 38) but does not experience pulsing (Parton et al., 2001). Equation (35) does quantify $r_{{\mathrm{NO}}_x/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$ as a function of Dr, but as a unitless ratio as expected.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{r}_{{\mathrm{NO}}_x/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}=\mathrm{15.2}\\ \text{(35)}& {\displaystyle }& {\displaystyle }+\left({\displaystyle \frac{\mathrm{35.5}\arctan \left(\mathrm{0.68}\mathit{\pi }\left(\left(\mathrm{10.0}\mathrm{Dr}\right)-\mathrm{1.86}\right)\right)}{\mathit{\pi }}}\right){\displaystyle }& {\displaystyle }{N}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}\left(\text{nitrification}{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\mathrm{g}\text{-}\mathrm{N}{\mathrm{h}}^{-\mathrm{1}}\right)=\\ \text{(36)}& {\displaystyle }& {\displaystyle }\mathrm{0.02}\left(R_{\mathrm{N}}\right)\left(\text{grid cell area}\right){\displaystyle }& {\displaystyle }{N}_{{\mathrm{NO}}_x}\left(\text{nitrification}{\mathrm{NO}}_x,\mathrm{g}\text{-}\mathrm{N}{\mathrm{h}}^{-\mathrm{1}}\right)=\\ \text{(37)}& {\displaystyle }& {\displaystyle }{r}_{{\mathrm{NO}}_x/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}\left(N_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}\right)P\left(l_{\mathrm{dry}}\right){\displaystyle }& {\displaystyle }{D}_{\mathrm{NO}}(\text{denitrification NO},\mathrm{g}\text{-}\mathrm{N}{\mathrm{h}}^{-\mathrm{1}})=\\ \text{(38)}& {\displaystyle }& {\displaystyle }{r}_{{\mathrm{NO}}_x/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}\left(D_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}\right)\end{array}$$

N₂O from denitrified NO₃ (R_D) is calculated using the partitioning function derived by Del Grosso et al. (2000) (Eq. 39). The ratio of N₂ to N₂O emitted as an intermediate during denitrification ($r_{{\mathrm{N}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$) is dependent on WFPS (Eq. 42) and the relative availability of NO₃ substrate and C for heterotrophic respiration (Eqs. 40–41). The C available for heterotrophic respiration in the surface soil layer (labile C) (Eq. 41) is taken from EPIC for agricultural biomes and from Xu et al. (2015) for nonagricultural biomes. f(NO₃:C) is controlled by variability in soil texture, accounted for by a factor k, which depends on soil diffusivity at field capacity as estimated in Del Grosso et al. (2000). Also, the NO₃ pool is updated at each time step when denitrification happens (Eq. 43). Equations (40)–(42) also quantify $r_{{\mathrm{N}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$ as a unitless ratio, while still accounting for the variables influencing these ratios.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{D}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}(\text{denitrification}{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\mathrm{g}\text{-}\mathrm{N}{\mathrm{h}}^{-\mathrm{1}})=\\ \text{(39)}& {\displaystyle }& {\displaystyle }\left({\displaystyle \frac{R_{\mathrm{D}}}{\mathrm{1.0}+r_{{\mathrm{N}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}}\right)\left(\text{grid cell area}\right)\text{(40)}& {\displaystyle }& {\displaystyle }{r}_{{\mathrm{N}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}=f\left({\mathrm{NO}}_{\mathrm{3}}:\mathrm{C}\right){\displaystyle }& {\displaystyle }f\left({\mathrm{NO}}_{\mathrm{3}}:\mathrm{C}\right)=\\ \text{(41)}& {\displaystyle }& {\displaystyle }\left\{\begin{array}{l}max\left(\mathrm{0.16}\left(k\right),\left(k\right)e^{-\mathrm{0.8}\left(\frac{{\mathrm{NO}}_{\mathrm{3}}}{\text{labile C}}\right)}\right)\text{if labile C}>\mathrm{0}\\ \mathrm{0.16}\left(k\right)\text{if labile C}\sim \mathrm{0}\end{array}\right.{\displaystyle }& {\displaystyle }f\left(\text{WFPS}\right)=\\ \text{(42)}& {\displaystyle }& {\displaystyle }max\left(\mathrm{0.1},\left(\mathrm{0.015}\left(\text{WFPS}\right(\text{as fraction})-\mathrm{0.32})\right)\right){\displaystyle }& {\displaystyle }{\mathrm{NO}}_{\mathrm{3}}\left(\mathrm{kg}\text{-}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}},\text{after denitrification}\right)=\\ \text{(43)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{R_{\mathrm{N}}}{K_{\mathrm{D}}}}+\left({\mathrm{NO}}_{\mathrm{3}}-{\displaystyle \frac{R_{\mathrm{N}}}{K_{\mathrm{D}}}}\right)\left(e^{-\left(K_{\mathrm{D}}\mathrm{d}t\right)}\right)\end{array}$$

HONO is emitted as an intermediate during nitrification and has been reported in terms of a ratio relative to NO for each of 17 ecosystems by Oswald et al. (2013). In the mechanistic scheme, the proportions of HONO relative to total NO_x for these 17 biomes were mapped to the closest 24 MODIS-type biome categories (Table A1) and then to the NLCD 40 types (HONO_f) with the mappings in Tables A2 and A3. This allows for consistency with sub-grid land use fractions from NLCD40. HONO emissions are further adjusted to reflect their dependence on WFPS (Oswald et al., 2013). The adjustment factor f_SWC reflects observations that HONO emissions rise linearly up to 10 % WFPS and then decrease until they are negligible around ∼40 % (Su et al., 2011; Oswald et al., 2013) (Eq. 45). Subsequently, total NO emission is a sum of nitrification NO emission, which is a difference of $N_{{\mathrm{NO}}_x}$ and S_HONO, and denitrification NO (Eq. 46). Similarly, total N₂O is a sum of $N_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$ (Eq. 36) and $D_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$ (Eq. 39). The canopy reduction factor (Sect. 2.1) is then applied to both S_HONO and S_NO (Eqs. 44 and 46). Finally, sub-grid-scale emission rates are aggregated for each grid cell.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{S}_{\mathrm{HONO}}=\left({\mathrm{HONO}}_f\right)\left(N_{{\mathrm{NO}}_x}\right)\\ \text{(44)}& {\displaystyle }& {\displaystyle }\left(f_{\mathrm{SWC}}\right)\mathrm{CRF}\left(\text{LAI, meteorology, biome}\right){\displaystyle }& {\displaystyle }{f}_{\mathrm{SWC}}\left(\text{soil water content adjustment factor to compute HONO}\right)=\\ \text{(45)}& {\displaystyle }& {\displaystyle }\left\{\begin{array}{l}{\displaystyle \frac{\left({\mathrm{HONO}}_f\right)\left(\text{WFPS}\right)}{\mathrm{0.1}}},\text{if}\left(\text{WFPS}\le \mathrm{0.10}\right)\\ \left(\text{assuming linear increase up to 10\hspace{0.17em}\% WFPS}\right)\\ {\displaystyle \frac{\left({\mathrm{HONO}}_f\right)\left(\mathrm{0.4}-\text{WFPS}\right)}{(\mathrm{0.4}-\mathrm{0.1})}},\text{if}\left(\text{WFPS}\le \mathrm{0.40}\right)\\ \mathrm{0},\text{if}\left(\text{WFPS}>\mathrm{0.40}\right)\end{array}\right.{\displaystyle }& {\displaystyle }{S}_{\mathrm{NO}}=\left\{\left(N_{{\mathrm{NO}}_x}-\left(\left({\mathrm{HONO}}_f\right)\left(N_{{\mathrm{NO}}_x}\right)\left(f_{\mathrm{SWC}}\right)\right)\right)+D_{\mathrm{NO}}\right\}\\ \text{(46)}& {\displaystyle }& {\displaystyle }\mathrm{CRF}\left(\text{LAI, meteorology, biome}\right)\end{array}$$

2.7 Model configurations

We obtained from the U.S. EPA a base case WRFv3.7-CMAQv5.1 simulation for 2011 with the settings and CONUS modeling domain described by Appel et al. (2017), who thoroughly evaluated its performance against observations. Here, we simulate only May and July to test the sensitivity of air pollution to soil N emissions during the beginning and middle of the growing season. Each episode is preceded by a 10-day spin-up period.

Table 2 summarizes the WRF-CMAQ modeling configurations settings. The simulations use the Pleim–Xiu Land Surface Model (PX-LSM) (Pleim and Xiu, 2003) and the Asymmetric Convective Mixing v2 (ACM2) planetary boundary layer (PBL) model. The modeling domain for CMAQ v5.1 covers the entire CONUS including portions of northern Mexico and southern Canada with 12 km resolution and a Lambert conformal projection. Vertically, we use 35 vertical layers of increasing thickness extending up to 50 hPa. Boundary conditions are provided by a 2011 global GEOS-Chem simulation (Bey et al., 2001).

Table 2 Modeling configuration used for the WRF-CMAQ simulations.

WRF simulations employed the same options as Appel et al. (2017) (summarized in Table 2). WRF outputs for meteorological conditions were converted to CMAQ inputs using MCIP version 4.2 (https://www.cmascenter.org, last access: 22 February 2019). Gridded speciated hourly model-ready emission inputs were generated using the Sparse Matrix Operator Kernel Emissions (SMOKE; https://www.cmascenter.org/smoke/, last access: 22 February 2019) version 3.5 program and the 2011 National Emissions Inventory v1. Biogenic emissions were processed in-line in CMAQ v5.1 using BEIS version 3.61 (Bash et al., 2016). All the simulations employed the bidirectional option for estimating the air–surface exchange of ammonia. We applied CMAQ with three sets of soil NO emissions: (a) standard YL soil NO scheme in BEIS; (b) updated BDSNP scheme for NO (Rasool et al., 2016) with new sub-grid biome classification; and (c) mechanistic soil N scheme for NO and HONO.

2.8 Observational data for model evaluation

To evaluate model performance for each of the three soil N cases, we employed regional and national networks: the EPA's Air Quality System (AQS; 2086 sites; https://www.epa.gov/aqs, last access: 22 February 2019) for hourly NO_x and O₃; the Interagency Monitoring of Protected Visual Environments (IMPROVE; 157 sites; http://vista.cira.colostate.edu/improve/, last access: 22 February 2019) and Chemical Speciation Network (CSN; 171 sites; https://www3.epa.gov/ttnamti1/speciepg.html, last access: 22 February 2019) for PM_2.5 nitrate (measured every third or sixth day); the Clean Air Status and Trends Network (CASTNET; 82 sites; http://www.epa.gov/castnet/, last access: 22 February 2019) for hourly O₃ and weekly aerosol PM species; and SEARCH network measurements (http://www.atmospheric-research.com/studies/SEARCH/index.html, last access: 22 February 2019) of NO_x concentrations in remote areas. NO₂ was also evaluated against tropospheric columns observed by the OMI aboard NASA's Aura satellite (Bucsela et al., 2013; Lamsal et al., 2014).

3 Results and discussion

3.1 Spatial distribution of soil NO, HONO, and N₂O emissions

Figure 3 compares the spatial distribution of soil N oxide emissions from the three schemes. The incorporation of EPIC fertilizer in BDSNP results in soil NO emission rates up to a factor of 1.5 higher than in YL, consistent with the findings of Rasool et al. (2016). Hudman et al. (2012) found nearly twice as large of a gap between BDSNP and YL in GEOS-Chem; the narrower gap here likely results from our use of sub-grid biome classification and EPIC fertilizer data (Rasool et al., 2016). The mechanistic scheme (Fig. 3c) generates emission estimates that are closer to the YL scheme but with greater spatial and temporal heterogeneity, reflecting its use of more dynamic soil N and C pools. The agricultural plains extending from Iowa to Texas with high fertilizer application rates have the highest biogenic NO and HONO emission rate, with obvious temporal variability between May and July (Fig. 3). In all of the schemes, soil N represents a substantial fraction of total NO_x emissions over many rural regions, especially in the western half of the country (Fig. S1 in the Supplement). However, the aggregated budget of soil NO is much less than anthropogenic NO_x from non-soil-related sources because fossil fuel use is concentrated in a limited number of urbanized and industrial locations. The percentage contribution of soil NO to total NO_x aggregated across the CONUS domain varied for May–July from 15 %–20 % for YL, 20 %–33 % for updated BDSNP, and 10 %–13 % for the mechanistic scheme.

Figure 3 Soil N oxide emissions on a monthly average basis for May (left) and July (right) 2011 for (a) the YL scheme (NO), (b) parameterized BDSNP scheme (NO), and (c) mechanistic scheme (NO + HONO).

Direct observations of soil emissions are sparse and most were reported decades ago. While the meteorological conditions will differ, these observations give us the best available indicator of the ranges of magnitudes of emission rates actually observed in the field. The sites encompass a variety of fertilized agricultural fields and fertilized and unfertilized grasslands (Bertram et al., 2005; Hutchinson and Brams, 1992; Parrish et al., 1987; Williams and Fehsenfeld, 1991, 1992; Martin et al., 1998). For fair comparison, the peak location or site was selected across a range of sites for a specific observation study and compared to the respective peak modeled value across sites or grids in the same spatial domain. Also, for comparison with natural unfertilized grassland observational studies based in Colorado, modeled estimates from nonagricultural grids only were selected. Overall, the YL scheme and the mechanistic scheme produce emissions estimates that are roughly consistent with the ranges of emission rates observed at each site (Table 3). By contrast, BDSNP tends to overestimate soil NO compared to these observations (Table 3).

Table 3 NO emission rates (ng N m⁻² s⁻¹) observed in field studies in agricultural and grassland locations, modeled by CMAQ with the three soil N schemes for May and July 2011. Observed and modeled values are from peak location or site within a range of values across sites.

^a Derived from SCIAMACHY NO₂ columns. ^b Mechanistic scheme estimates are NO + HONO emission rates.

Table 3 also shows opposing trends for May and July soil NO estimates between YL or BDSNP and the mechanistic scheme for Iowa and South Dakota fertilized fields that make up a significant part of the Corn Belt in the US. For these regions, soil NO tends to be higher in July than in May in YL and BDSNP, but lower in July in the mechanistic scheme (Table 3). The US Corn Belt has the most synthetic N fertilizer application in April (Wade et al., 2015), which can explain the high soil NO emissions in May that decline in July. N₂O emissions have been particularly observed to be highest during May–June after April N fertilizer application in the US Corn Belt, with a decline thereafter (Griffis et al., 2017). This is further confirmed in our estimates for soil N₂O emissions from the mechanistic scheme, for which May estimates are higher than in July and the maximum emissions are observed in the Iowa Corn Belt (Fig. 4). However, unlike NO_x emissions, for N₂O no background conditions or emission inventories are in place in CMAQ's chemical transport model, so comparisons with ambient observations are not yet possible.

Figure 4 Soil N₂O emissions on a monthly average basis for May (a) and July (b) 2011 estimated from the mechanistic scheme.

3.2 Evaluation with PM_2.5, ozone, and NO_x observations

Model results with the three soil N schemes are compared with observational data from IMPROVE and CSN monitors for the PM_2.5 NO₃ component, AQS monitors for NO_x and ozone, and CASTNET monitors for ozone. Both YL and the new mechanistic scheme exhibit similar ranges of bias for these pollutants (see Figs. S2, S3, S4, S5, and S6 in the Supplement). Use of the mechanistic scheme in place of YL changes soil N emissions by less than 25 ng N m⁻² s⁻¹ in most regions, corresponding to NO_x concentration changes of less than 1 ppb (Fig. 5). CASTNET and IMPROVE monitors tend to be more remote than AQS and CSN monitors, many of which are located in urban regions.

Figure 5 Total NO_x (NO + NO₂) concentration sensitivity (right) to changes in soil NO_x emissions (left) on a monthly average basis for May (a, c) and July (b, d) 2011 when switching from the YL scheme (NO) to the mechanistic scheme (NO + HONO).

At AQS monitors, switching between soil N schemes changes MB for O₃ by up to ∼1.5 ppb (Fig. 6), whereas the absolute MB of models versus observations is up to ∼10 ppb (Fig. S2). For NO_x, the maximum difference in MB between soil N schemes is ∼0.4 ppb (Fig. 7) compared to a maximum absolute MB of ∼10 ppb between model and observations (Fig. S3). For CASTNET monitors, the differences in MB for O₃ between soil N schemes can reach a maximum of ∼1.5 ppb (Fig. 8) compared to the 6 ppb maximum absolute MB of models versus observations (Fig. S4). Similarly, for IMPROVE PM_2.5 NO₃, the maximum difference in MB between soil N schemes is ∼0.06 µg m³ (Fig. 9) compared to the maximum absolute MB of 0.4 µg m³ (Fig. S5). For CSN PM_2.5 NO₃, the maximum MB difference between soil N schemes is ∼0.1 µg m³ (Fig. 10) compared to the maximum absolute MB of ∼50 µg m³ (Fig. S6). Similar trends are observed for both May and July as illustrated in Figs. 6–10.

Overall, the mechanistic scheme tends to reduce CMAQ's positive biases for pollutants across the Midwest and eastern US, whereas BDSNP worsens overestimations in these regions for both May and July 2011 (Figs. 6–10). In addition, the negative bias in the difference means less bias compared to observations (Figs. 6–10). One reason for the differences is that the mechanistic scheme recognizes dry conditions in unirrigated fields in these regions, whereas the low WFPS threshold in BDSNP (θ=0.175 (m³ m⁻³)) treats most of these regions as wet and thus higher emitting.

Figure 6 Change in average monthly mean bias (MB) of the Community Multiscale Air Quality (CMAQ) model evaluated against the EPA Air Quality System (AQS) O₃ observations for May (a, b) and July (c, d) 2011 when switching to the mechanistic (a, c) or BDSNP (b, d) scheme from YL.

Figure 7 Change in average monthly MB of CMAQ evaluated against EPA AQS NO_x observations for May (a, b) and July (c, d) 2011 when switching to the mechanistic (a, c) or BDSNP (b, d) scheme from YL.

Figure 8 Change in average monthly MB of CMAQ evaluated against the EPA Clean Air Status and Trends Network (CASTNET) O₃ observations for May (a, b) and July (c, d) 2011 when switching to the mechanistic (a, c) or BDSNP (b, d) scheme from YL.

Figure 9 Change in average monthly MB of CMAQ evaluated against Interagency Monitoring of Protected Visual Environments (IMPROVE) PM_2.5 NO₃ observations for May (a, b) and July (c, d) 2011 when switching to the mechanistic (a, c) or BDSNP (b, d) scheme from YL.

Figure 10 Change in average monthly MB of CMAQ evaluated against Chemical Speciation Network (CSN) PM_2.5 NO₃ observations for May (a, b) and July (c, d) 2011 when switching to the mechanistic (a, c) or BDSNP (b, d) scheme from YL.

3.2.1 Evaluation with South Eastern Aerosol Research and CHaracterization (SEARCH) network NO_x measurements

We analyzed how the choice of soil NO parameterization affects NO_x concentrations in nonagricultural regions by using SEARCH network measurements (http://www.atmospheric-research.com/studies/SEARCH/index.html, last access: 22 February 2019). Six SEARCH sites located in the southeastern US are evaluated for May and July 2011: Gulfport, Mississippi (GFP), an urban coastal site ∼1.5 km from the shoreline; Pensacola, an outlying (aircraft) landing field (OLF) remote coastal site near the gulf ∼20 km inland; Atlanta, Georgia (Jefferson Street, JST), and North Birmingham, Alabama (BHM) – both urban inland sites; and Yorkville, Georgia (YRK), and Centreville, Alabama (CTR) – remote inland forest sites.

Figure 11 Comparison of average monthly (May and July 2011) MB for CMAQ NO_x with (a) YL, (b) BDSNP parameterized, and (c) mechanistic schemes compared to South Eastern Aerosol Research and CHaracterization (SEARCH) NO_x observations in nonagricultural remote regions.

Across the southeastern US during these episodes, BDSNP estimated higher emissions than YL and the mechanistic scheme estimated lower emissions (Fig. 3). Also, CMAQ with each scheme overestimated NO_x observed at each SEARCH site (Fig. 11). Thus, shifting from YL to BDSNP worsens mean bias (MB) for NO_x, while the mechanistic scheme reduces MB. The impacts are most pronounced at the rural Centreville site (Fig. 11).

3.3 Evaluation with OMI satellite NO₂ column observations

Tropospheric NO₂ columns observed by OMI and available publicly at the NASA archive (http://disc.sci.gsfc.nasa.gov/Aura/data-holdings/OMI/omno2_v003.shtml, last access: 22 February 2019; Bucsela et al., 2013; Lamsal et al., 2014) are used to evaluate the performance of CMAQ under the three soil NO_x schemes. To enable a fair comparison, the quality-assured and quality-checked (QA ∕ QC) clear-sky (cloud radiance fraction <0.5) OMI NO₂ data are gridded and projected to our CONUS domain using ArcGIS 10.3.1. CMAQ NO₂ column densities in molecules per cm² are generated from CMAQ through vertical integration using the variable layer heights and air mass densities in these tropospheric layers. These NO₂ column densities are then extracted for 13:00–14:00 local time across the CONUS domain to match the time of OMI overpass measurements.

We compared CMAQ-simulated tropospheric NO₂ columns with OMI data for four broad regions that showed the highest sensitivity to the soil N schemes. For May 2011, the mechanistic scheme produces higher estimates of NO₂ than YL in the western US and Texas, with lower estimates in the rest of the agricultural Great Plains. In July, however, the mechanistic scheme produces lower estimates than YL in each of these regions, but the differences are narrower than in May (Fig. 12). Switching from YL to our updated mechanistic scheme improved agreement with OMI NO₂ columns in the western US (for May only), Montana, North and South Dakota, North and South Carolina and Georgia (July only), and Oklahoma and Texas (red boundaries). However, switching from YL to the mechanistic scheme worsens underpredictions of column NO₂ in the rest of the Midwest (black boundaries) during both May and July (Figs. 12 and 13). The mechanistic scheme improves model performance in the southeastern US and many portions of the central and western US (Table 4). Overestimation is exhibited for the eastern US across all soil N schemes and can be attributed more to the current emission inventory in CMAQ overestimating NO₂ vertical column density in this region of CONUS (Kim et al., 2016). For Texas and Oklahoma, the mechanistic scheme performs better than YL but still underestimates OMI observations in May and performs well in July (Fig. 13).

Figure 12 Impact of switching from the YL scheme to the mechanistic scheme on CMAQ tropospheric NO₂ column density at NASA's Ozone Monitoring Instrument (OMI) overpass time (13:00–14:00 local time) on a monthly average (May and July 2011) basis.

Table 4 Statistical performance of the CMAQ modeled (with YL, updated BDSNP, and mechanistic schemes) tropospheric NO₂ column for May 2011 with OMI NO₂ observations for sensitive sub-domains for CONUS.

Underestimates of soil N in some regions with an abundance of animal farms, such as parts of Colorado, New Mexico, north Texas, California, the northeast US, and the Midwest, may be attributed to the lack of representation of farm-level manure N management practices, in which manure application can exceed the EPIC estimate of optimal crop demand. Farms in the vicinity of concentrated animal units often apply N in excess of the crop N requirements as part of the manure management strategy, typically increasing the N emissions (Montes et al., 2013). The USDA has reported that confined animal units or livestock production correlates with increasing amounts of farm-level excess N (Kellogg et al., 2000; Ribaudo et al., 2016). Model representations of these practices are needed to better estimate the impact of nitrogen in the environment.

4 Conclusions

Our implementation of a mechanistic scheme for soil N emissions in CMAQ provides a more physically based representation of soil N than previous parametric schemes. To our knowledge, this is the first time that soil biogeochemical processes and emissions across a full range of nitrogen compounds have been simulated in a physically realistic manner in a regional photochemical model. Our mechanistic scheme directly simulates nitrification and denitrification processes, allowing it to consistently estimate soil emissions of NO, HONO, NH₃, and N₂O (Figs. 1 and 2). The mechanistic scheme also updates the representation of the dependency of soil N on WFPS by utilizing parameters like water content at saturation, wilting point, and field capacity and their impact on gas diffusivity (Del Grosso et al., 2000; Parton et al., 2001).

Overall, the magnitudes of soil NO_x emissions predicted by the mechanistic scheme are similar to those predicted by the YL parametric scheme and smaller than those predicted by the BDSNP scheme. In dry conditions, soil NO has been shown to be the highest compared to wet conditions with the lowest, explained by sustained high nitrification rates due to high gas diffusivity in dry conditions (Homyak and Sickman, 2014). Arid soils or dry seasons with adequate soil N due to asynchrony between soil C mineralization and nitrification have been shown to shut down plant N uptake through high gas diffusivity, causing NO emissions to increase (Evans and Burke, 2013; Homyak et al., 2016). The mechanistic scheme exhibits this spatial variability in soil NO depending on dry or wet conditions, since it accounts for their dependence on soil moisture and gas diffusivity, as well as the C and N cycling that leads to adequate soil N.

Spatial patterns of NO_x emissions differ across the schemes and episodes (Fig. 3), but generally show the highest emissions in fertilized agricultural regions. During the episodes considered here, Texas experienced severe to extreme drought, while parts of the northeast and Pacific Northwest were unusually wet (http://www.cpc.ncep.noaa.gov/products/analysis_monitoring/regional_monitoring/palmer/2011/, last access: 22 February 2019). Testing for other time periods is needed to see how results differ during different seasons and as drought conditions vary. Model evaluation will also depend on the meteorological model's skill in capturing dry and wet conditions.

Figure 13 Comparison of average monthly (May and July 2011) OMI NO₂ column densities with CMAQ tropospheric NO₂ column density using YL, BDSNP, and mechanistic schemes. Regions are depicted in Fig. 12.

The lower emissions of the mechanistic scheme reduce the overprediction biases for ground-based observations of ozone and PM nitrate that had been reported by Rasool et al. (2016) for the BDSNP scheme (Figs. 6–10). The mechanistic scheme reduced overpredictions of NO_x concentrations at SEARCH sites in the southeastern US (Fig. 11). However, changes in performance for simulating satellite observations of NO₂ columns were mixed (Figs. 12–13). The underestimation of NO₂ by CMAQ with the mechanistic scheme in agricultural regions of the Midwest may be partially attributed to neglecting manure management practices from livestock operations. In the US, 60 % of nitrogen from manure produced on animal feedlot operations cannot be applied back to the same land because it is in “excess” of USDA advised agronomic rates. Most US counties with animal farms have adequate crop acres not associated with animal operations, but these are within the county, on which it is feasible to spread the excess manure at agronomic rates at certain additional cost. However, 20 % of the total US on-farm excess manure nitrogen is produced in counties with insufficient cropland for its application at agronomic rates (Gollehon et al., 2001). For areas without adequate land, alternatives to local land application such as energy production (for example, biofuel) are needed. In the absence of such a mitigation strategy, excess manure N applied on soil contributes to reactive N emissions and leaching (Ribaudo et al., 2003, 2012).

Although this work represents the most process-based representation of soil N ever introduced to a regional photochemical model, limitations remain. EPIC still lacks a complete representation of farming management practices like excess N applied as part of a nutrient management strategy for livestock, which can increase soil N pools and associated emissions. Developing and evaluating these models to address management decisions is challenging as they are often regionally specific and based on expert knowledge including regional and global economics and biogeochemical processes that have yet to be codified into a predictive system. Some aspects of soil N biogeochemistry remain insufficiently understood, especially as they relate to HONO emissions. Nevertheless, the mechanistic approach introduced here will make it possible to incorporate future advancements in understanding C and N cycling processes.

For future work, there is a need for more accurate representation of actual farming practices beyond the generalizations made by the EPIC model. Model development should be continued to better constrain N sources such as rock weathering, which are still ignored for estimating soil N emissions. Recently, Houlton et al. (2018) postulated that bedrock weathering can contribute an additional 6 %–17 % to global inorganic soil N for different natural biomes. There is also a need for more field observations of soil N emissions to better evaluate the spatial and temporal patterns simulated by the models.

Appendix A

Table A1 List of 24 MODIS soil-biome-based Cmic, Nmic, and HONO_f emission factors (%) derived from Xu et al. (2013) and Oswald et al. (2013).

^* A – equatorial, B – arid, C – warm temperature, D – snow, E – polar

Table A2 Mapping table to create the MODIS 24 soil biome map based on NLCD40 MODIS land cover categories for updated BDSNP parameterization.

^* NLCD categories 18 and 19 were mapped as MODIS category 1 (water) in Rasool et al. (2016), which have been corrected here.

Table A3 Microbial ∕ organic biomass C and N % and $\mathrm{HONO}/{\mathrm{N}}_{{\mathrm{NO}}_x}$ % mapped to the respective NLCD40 MODIS land cover categories based on Xu et al. (2013) estimates.

^a NLCD classes 26 and 27 consisting of mostly rocks. ^b Cmic and Nmic for US croplands classified under NLCD classes 37 and 38 are kept as zero to prevent double counting, as they are accounted for by EPIC N data.