Articles | Volume 23, issue 17
Research article
31 Aug 2023
Research article |  | 31 Aug 2023

Direct observations of NOx emissions over the San Joaquin Valley using airborne flux measurements during RECAP-CA 2021 field campaign

Qindan Zhu, Bryan Place, Eva Y. Pfannerstill, Sha Tong, Huanxin Zhang, Jun Wang, Clara M. Nussbaumer, Paul Wooldridge, Benjamin C. Schulze, Caleb Arata, Anthony Bucholtz, John H. Seinfeld, Allen H. Goldstein, and Ronald C. Cohen

Nitrogen oxides (NOx) are principle components of air pollution and serve as important ozone precursors. As the San Joaquin Valley (SJV) experiences some of the worst air quality in the United States, reducing NOx emissions is a pressing need, yet quantifying current emissions is complicated due to a mixture of mobile and agriculture sources. We performed airborne eddy covariance flux measurements during the Re-Evaluating the Chemistry of Air Pollutants in California (RECAP-CA) field campaign in June 2021. Combining footprint calculations and land cover statistics, we disaggregate the observed fluxes into component fluxes characterized by three different land cover types. On average, we find emissions of 0.95 mg N m−2 h−1 over highways, 0.43 mg N m−2 h−1 over urban areas, and 0.30 mg N m−2 h−1 over croplands. The calculated NOx emissions using flux observations are utilized to evaluate anthropogenic emissions inventories and soil NOx emissions schemes. We show that two anthropogenic inventories for mobile sources, EMFAC (EMission FACtors) and FIVE (Fuel-based Inventory for Vehicle Emissions), yield strong agreement with emissions derived from measured fluxes over urban regions. Three soil NOx schemes, including the MEGAN v3 (Model of Emissions of Gases and Aerosols from Nature), BEIS v3.14 (Biogenic Emission Inventory System), and BDISNP (Berkeley–Dalhousie–Iowa Soil NO Parameterization), show substantial underestimates over the study domain. Compared to the cultivated soil NOx emissions derived from measured fluxes, MEGAN and BEIS are lower by more than 1 order of magnitude, and BDISNP is lower by a factor of 2.2. Despite the low bias, observed soil NOx emissions and BDISNP present a similar spatial pattern and temperature dependence. We conclude that soil NOx is a key feature of the NOx emissions in the SJV and that a biogeochemical-process-based model of these emissions is needed to simulate emissions for modeling air quality in the region.

1 Introduction

Nitrogen oxides (NOxNO+NO2) are important trace gases that affect both the gas and aerosol phases of tropospheric chemistry. NOx regulates the concentrations of the primary atmospheric oxidant, hydroxyl radicals (OH), and serves as the catalyst for the formation of ozone (O3). NOx also affects the formation of inorganic nitrate aerosol through the production of nitric acid (HNO3) and organic nitrates (RONO2) and plays a role in secondary organic aerosol (SOA) production. NOx, O3, and aerosol are all detrimental to human health, triggering respiratory diseases (Kampa and Castanas2008; Hakeem et al.2016) and leading to premature death (Lelieveld et al.2015).

NOx is predominantly emitted from anthropogenic sources, including light- and heavy-duty transportation, fuel combustion, and biomass burning. Among these sectors, transportation is the largest in the United States (EPA2016). Strict regulations have been implemented to control NOx emissions. Three-way catalysts have effectively reduced emissions from gasoline-powered passenger vehicles. The application of emissions control systems on coal power plants has reduced NOx emissions (De Gouw et al.2014). The California Air Resources Board (CARB) has proposed the Heavy-Duty Engine and Vehicle Omnibus Regulation and associated amendments and a target for a 90 % reduction in per-vehicle heavy-duty NOx emissions by 2031 (CARB2016). The regulation of mobile sources leads to an increasing importance of natural NOx sources, such as lightning and soil emissions. Soil NOx is released as a byproduct of microbial nitrification and denitrification (Andreae and Schimel1990). While the biogeochemistry of soil NOx emissions is well established, this biogenic source involves a complex interaction of soil microbial activity, namely the soil nitrogen (N) content. Besides, agriculture activities, such as the use of fertilizers, lead to a substantial enhancement of soil NOx emissions (Phoenix et al.2006).

Currently, the San Joaquin Valley (SJV) in California experiences some of the most severe air pollution in the United States. The SJV cities, Visalia, Fresno, and Bakersfield are among the top 10 most polluted cities for both ozone and particulate matter (American Lung Association2020). In order to implement appropriate emissions control efforts, identifying the contribution of different NOx emissions is particularly important for the SJV, as it features a complex mixture of emissions from fuel combustion and soil emissions associated with agriculture. The contribution of soil NOx emissions remains highly uncertain. While Guo et al. (2020) attribute approximately 1.1 % of anthropogenic NOx emissions in California to soil NOx, Almaraz et al. (2018) argued that due to growing N fertilizer use, the SJV has soil NOx emissions of 24 kg N ha−1 yr−1, contributing to 20 %–51 % of the NOx budget of the entire state of California. Similarly, Sha et al. (2021) estimated that 40.1 % of the total NOx emissions over California in July 2018 are from soils.

Airborne eddy covariance (EC) flux measurements provide a powerful tool to investigate the emissions strength of atmospheric constituents at landscape scales. It has been applied to assess the surface exchanges of greenhouse house gases (GHGs), including CO2 and methane (CH4) (Mauder et al.2007; Yuan et al.2015; Sayres et al.2017; Hannun et al.2020). In recent years, it has been extended to the study of emissions of volatile organic compounds and NOx over a megacity (Karl et al.2009; Vaughan et al.2021), vegetation (Karl et al.2013; Misztal et al.2014; Wolfe et al.2015; Kaser et al.2015; Yu et al.2017; Gu et al.2017), and shale gas production regions (Yuan et al.2015). Compared to the traditional EC measurements from instruments mounted at a fixed location on a tower, wavelet-based airborne EC measurements allow for larger spatial assessment and are well suited to regions with inhomogeneous and non-stationary source distributions (Sühring et al.2019).

In this study, we present airborne EC flux measurements obtained during seven flights of a Twin Otter aircraft over the San Joaquin Valley in California. Companion studies of NOx emissions over Los Angeles (Nussbaumer et al.2023) and volatile organic compound (VOC; Pfannerstill et al.2023) and GHG fluxes (Schulze et al.2023) will be presented separately. We utilize continuous wavelet transformation to calculate the NOx flux (Sect. 3). In conjunction with the footprint calculations and land classification, we explore the spatial heterogeneity of NOx emissions and identify component fluxes from the highway, urban, and soil land types (Sect. 4). We also utilize the NOx emissions derived from flux measurements to evaluate anthropogenic emissions inventories and soil NOx schemes (Sect. 5).

2 Measurements

The airborne EC flux measurements were conducted on a Twin Otter research aircraft operated by the Naval Postgraduate School (NPS) during the Re-Evaluating the Chemistry of Air Pollutants in California (RECAP-CA) field campaign. The RECAP-CA field campaign was conducted between 1 and 22 June in California, including 7 d of measurements over the San Joaquin Valley and 9 d of measurements over Los Angeles. The flight path was designed with long straight legs to ensure the good quality of the flux measurements (Fig. S1 in the Supplement; Karl et al.2013). The aircraft flew slowly, at an airspeed of 50–60 m s−1, and cruised at a low height of  300 m above ground. The aircraft took off at  11:00 local time (LT) at Hollywood Burbank Airport and landed at  18:00 LT.

The standard instruments aboard the aircraft are described in Karl et al. (2013) and include total and dew point temperature, barometric and dynamic pressures, wind direction and wind speed, total airspeed, slip and attack angles, GPS latitude, GPS longitude, GPS altitude, pitch, roll, and heading. These measurements are at 10 Hz temporal resolution. VOCs were measured at 10 Hz time resolution by Vocus proton transfer reaction time-of-flight mass spectrometer (Vocus PTR-ToF-MS), as described in Pfannerstill et al. (2023). Mixing ratios of NOx were measured at 5 Hz frequency, using a custom-built three-channel thermal dissociation–laser-induced fluorescence (TD–LIF) instrument. The multipass LIF cells, fluorescence collection, long-pass wavelength filtering (for λ>700 nm), and photon counting details have been previously described (Thornton et al.2000; Day et al.2002; Wooldridge et al.2010). Details specific to this implementation are described below.

Air was sampled from the aircraft community inlet through PFA (perfluoroalkoxy) Teflon tubing at a rate of ∼6 L min−1 and split equally between the three instrument channels. Each channel measured NO2 by laser-induced fluorescence, utilizing a compact green laser (Spectra-Physics Explorer One XP; 532 nm). The laser was pulsed at 80 kHz, and the 1.7 W average power was split between the three cells. Earlier versions of the instrument used a dye laser that was tuned on and off a narrow rovibronic NO2 resonance at 585.1 nm. Experience over a wide variety of conditions had demonstrated that the offline signal did not depend on the sample, other than from aerosol particles, and that could be eliminated by adding a Teflon membrane filter. Moving to nonresonant excitation at 532 nm provided full-time coverage at 5 Hz, along with lower complexity and more robust performance of the laser system. Maintaining the LIF cells at low pressure (∼0.4 kPa) was no longer required to avoid line broadening but was still desirable to extend the NO2 fluorescence lifetime for time-gated photon counting in order to reject prompt laser scatter. Instrument zeros were run using ambient air scrubbed of NOx every 20 min in flight to correct for any background drift during the flights. In addition, calibrations were performed in-flight every 60 min, using a NO2 in N2 calibration cylinder (Praxair; 5.5 ppm; Certified Standard grade) diluted with scrubbed air.

NO2 was measured directly in the first channel, with the sample passing only through a particle filter and a flow-limiting orifice before the cell. NOx was measured in the second by adding O3 (generated with 184.5 nm light and a flow of scrubbed and dried air) to convert NO to NO2 before detection. A 122 cm length of 0.4 cm internal diameter tubing served as the O3+NO reactor, providing 4 s of reaction time before the orifice. The third channel was used to measure the sum of higher nitrogen oxides (e.g., organic nitrates and nitric acid) by thermal dissociation to NO2 with an inline oven (∼500C) before LIF detection.

3 Flux and footprint calculation

3.1 Pre-processing

The observed 10 Hz vertical wind speeds are downscaled to 5 Hz in order to match the time resolution of the NOx measurements. The full observation data set breaks into segments with continuous wind and NOx measurements. The segment window is selected if the length is larger than 10 km and if the height variation is less than 200 m. We also filter out measurements when aircraft roll angles are larger than 8 to avoid perturbation in the vertical wind due to aircraft activity. While most of the measurements are within the planetary boundary layer (PBL), the airplane occasionally arose above the boundary layer, and these observations above PBL are removed in later analysis. The PBL heights are determined using the sharp gradient in the dew point, water concentration, toluene concentration, and temperature at the soundings conducted during the voyage, and we interpolate the PBL heights to the full duration of the flight. The PBL heights agree well against the hourly PBL heights from the High-Resolution Rapid Refresh (HRRR) product (Fig. S2).

We adjust for the lag time between the meteorology measurements and the TD–LIF measurements by shifting the time of the TD–LIF observation within the time window of ±4 s, until the covariance with the vertical wind speed is maximized (Fig. S4). As the time lag is assumed to be due to differences in the clocks of the two instruments and the transit time of air through the TD–LIF instrument, we assume that the lag time for each flight is constant. We use the median lag time from each flight for all segments collected on the same day.

3.2 Continuous wavelet transformation

The continuous wavelet transformation (CWT) parameterization decomposes the time series (x(t)) into a range of frequencies and represents it as the convolution of the time series with a wavelet function (Torrence and Compo1998).


where W(a,b) is the wavelet coefficient. ψa,b*(t) is the wavelet function, which is based on a “parent” wavelet ψ0, and is adjusted with a transition parameter b and a scale parameter a. The transition parameter determines the location of the “parent” wavelet, and the scale parameter defines the frequency. We use the Morlet wavelet as the “parent” wavelet.

(3) ψ 0 = π - 1 / 4 e 6 i η e - η 2 / 2

The Morlet wavelet has been widely applied to represent turbulence in the atmosphere due to a reasonable localization in the frequency domain and a good ability with respect to edge detection (Schaller et al.2017).

Time domain scales are increased linearly with the increment of the time resolution (δt, 0.2 s), and N is the number of data points. Frequency domain scales are represented by an exponential array of scale parameters aj with the increment δj of 0.25 s. The smallest frequency scale is the Nyquist frequency, which is twice the time resolution (0.4 s).


For two simultaneous time series of NOx (Wc(a,b)) and vertical wind speed (Ww(a,b)), we obtain the wavelet cross-spectrum, following Eq. (6). The Morlet-wavelet-specific reconstruction factor Cδ is 0.776. We then sum up over the full frequency scales to yield a time series of the flux (Eq. 7).


Figure 1 exhibits an example of CWT flux calculation. Figure 1a shows the normalized NOx and vertical wind speed in a straight segment of ∼50 km. The normalization is realized by subtracting out the average followed by dividing the standard deviation of a scalar time series. Both time series are decomposed using a CWT algorithm to yield the cross-power spectrum shown in Fig. 1b. Due to the finite length in time, the wavelet power spectrum is prone to higher uncertainties closer to the edge (Mauder et al.2007). The regions of the wavelet power spectrum in which the edge effects are the largest are identified as the cone of influence (COI). Data points containing >80 % spectral power within the cone of influence are removed for quality control. The power spectrum is then integrated over all frequencies to the time series of the NOx flux (Fig. 1c). We processed the integrated fluxes as 2 km moving averages to address the influence of large-scale turbulence and then resampled them at 500 m.

Figure 1(a) The variance in the normalized NOx and vertical wind speed. (b) Frequency and time-resolved wavelet power spectrum, with the cone of influence shown as a dotted black line. (c) The integrated fluxes from the raw data points are shown in black, and the fluxes after the moving average and COI filtering have been applied are shown in green. (d) The distribution map of flux resampled at 500 m. The black lines show the 90th percentiles of the footprints, and the thick black line denotes the contours of all footprints. Background map credits: Esri, HERE, Garmin, NGA, USGS, and NPS.

3.3 Footprint calculation

The footprint describes the contribution of surface regions to the observed airborne flux. We use the KL04-2D parameterization (Kljun et al.2004) to calculate a space-resolved footprint map. This KL04-2D parameterization is developed from a 1-D backward Lagrangian stochastic particle dispersion model (Kljun et al.2004). Metzger et al. (2012) implemented a Gaussian cross-wind distribution function to resolve the dispersion perpendicular to the main wind direction. The input parameters include the height of the measurements, standard deviation of the horizontal and vertical wind speed, horizontal wind direction, boundary layer height, surface roughness length, and friction velocity. We obtain the surface roughness length and friction velocity from the HRRR product. For each flux observation, we calculate the footprint map at the spatial resolution of 500 m and then extract the 90 % contour. Figure 1d depicts the 90 % KL04-2D footprint contours of observations resampled to 500 m in one segment. Each footprint contour is aligned with the horizontal wind direction and is transformed into a geographic coordinate space.

3.4 Filter out NOx fluxes impacted by the off-road vehicle emissions

It is worth noting that croplands include not only soil NOx emissions but also the off-road vehicle emissions. Erroneously attributing the NOx from off-road vehicle emissions to soil NOx emissions leads to a high bias. While trimethylbenzene was observed during RECAP-CA field campaign, Pfannerstill et al. (2023) presented the trimethylbenzene fluxes using the same algorithm described in Sect. 3.2. The trimethylbenzene fluxes are interpolated to match the NOx fluxes in time and are utilized as an indicator of off-road vehicle emissions over croplands (Tsai et al.2014). The trimethylbenzene fluxes are categorized into two groups; the first group presents footprints covering croplands exclusively, and the second group presents footprints with mixed land cover types. Shown in Fig. S9, the trimethylbenzene flux is much lower over croplands, with a median of 0.003 mg m−2 h−1 compared to a median of 0.009 mg m−2 h−1 over mixed land cover types, which includes highway and urban areas. No correlation between trimethylbenzene flux and NOx flux is found over croplands. Among all observations over croplands, we identify those with trimethylbenzene fluxes larger than 0.02 mg m−2 h−1, which consist of 7 % of the total data points and are impacted by the off-road vehicle emissions, and then filter out them in the later analysis. We also vary the threshold of the trimethylbenzene flux between 0.005 and 0.04 mg m−2 h−1 and conclude that the choice of the threshold does not influence the results.

3.5 Vertical divergence

Extrapolating the airborne flux to surface flux should account for the vertical divergence. The vertical divergence is a result of multiple processes, including net in situ production or loss, storage, and horizontal advection.

To investigate the impact of vertical divergence, the flight route includes three vertically stacked racetracks during which the segments are close to each other in space but vary in height. After removing the legs that fail the quality control, only one racetrack measurement carried out between 14:20–15:10 LT on 8 June presented qualified flux segments, and the vertical distribution of fluxes is shown in Fig. S4. No consistent increase or decrease in the fluxes with increasing height is detected during the racetrack in this study because the vertical divergence is hampered by emissions heterogeneity. As shown in Fig. S7, the footprint map for each segment at various altitudes covers regions with high heterogeneity. Therefore, we use an alternative approach for calculating the vertical divergence. Instead of extracting racetrack measurements, we collect a subset of flux measurements during the whole field campaign based on the footprint coverage. Only fluxes with footprints covering croplands exclusively are included to avoid emissions heterogeneity. As shown in Fig. 2, we calculate the ratio of measurement heights relative to the PBL height (z/zi), and 98 % of selected fluxes are located within 70 % of the PBL height and are divided into seven bins of z/zi with uniform width. We then perform a linear fit for the binned median fluxes versus z/zi to calculate the vertical correction factor C=slopeintercept. This correction factor is used to linearly extrapolated the fluxes at the measurement height (Fz) to fluxes at the surface (F0) (Eq. 8). After the vertical divergence correction, the surface fluxes are on average 26 % higher than the fluxes at the measurement heights.

(8) F 0 = F x 1 + C z z i

Figure 2Vertical profiles of measured fluxes above croplands during RECAP-CA field campaign binned by the ratio of measurement height and PBL height (z/zi). The points represent the median flux within each bin, and the error bars represent the standard deviation. The dashed red line shows a linear fit for median fluxes versus relative height.


3.6 Data qualify control and uncertainty analysis

The flux detection limit does not only depend on the signal-to-noise ratio of the NOx measurement but also varies with wind speed and atmospheric stability. Following Langford et al. (2015), we calculate the limit of detection (LoD) of flux before the moving average and spatial average are applied. For each segment, the observed NOx is replaced with a white noise time series and is then fed into the CWT to yield the corresponding time series of the noisy flux. The random error affecting the flux (σNOx,noise) is defined as the standard deviation of this noise-derived flux, and LoD is defined as 2×σNOx,noise (95th confidence level). Among 142 segments, Fig. 3a shows the distribution of flux LoD among 142 segments. The LoDs range from 0.02 to 0.30 mg N m−2 h−1, and the average LoD is 0.10 mg N m−2 h−1. To obtain a better constraint on the flux quality, we compare the LoD against the time series of flux in each segment and filter out 18 segments in which the whole time series is below the LoD. The flux calculation using CWT introduces uncertainty from a variety of sources. We describe the systematic errors and random errors, following the work of Wolfe et al. (2018). Systematic errors arise from the undersampling of high-frequency and low-frequency ranges. The CWT algorithm fails to resolve a frequency higher than the Nyquist frequency. Due to the high temporal resolution (5 Hz), we expect a minimal loss at the high-frequency limit (Fig. S5). The upper limit of the systematic error associated with a low frequency is calculated using Eq. (9) (Lenschow et al.1994).

(9) SE 2.2 z z i 0.5 z i L

z and L are the measurement heights and the length of segments, respectively. zi values are the boundary layer heights from HRRR. We calculate the low-frequency error ranges from 1 % to 5 %.

Figure 3(a) The distribution of the segment-based NOx flux detection limit (LoD). (b) The distribution of total uncertainty in the NOx flux at 500 m spatial average.


Random errors arise from the noise in the instrument (REnoise) and the noise in turbulence sampling (REturb), which are calculated using Eqs. (10) and (11) (Wolfe et al.2018; Lenschow et al.1994).

(10) RE noise = σ NO x , noise 2 σ w 2 N

(11) RE turb F 1.75 z z i 0.25 z i L 0.5

z, L, and zi values are the same as Eq. (9), and σw2 is the variance of the vertical wind speed. Note that REnoise assumes that the noise in each time step is uncorrelated; therefore, we ignore the moving average step in the uncertainty calculation, and N denotes the number of points used to yield each 500 m spatially averaged flux.

Utilizing a constant lag time introduces an additional source of uncertainty. We estimate the uncertainty by comparing the calculated fluxes using segment-specific and constant lag times across all segments for which specific lag times are available. As shown in Fig. S4, the difference is less than 25 % for 90 % of the data. Therefore, we attribute an uncertainty of 25 % due to the lag time correction (RElag). While we believe that this error is unphysical and that a single lag time is more appropriate, we include it to be conservative in our estimate of the uncertainties.

Estimating the uncertainty caused by the correction of vertical divergence is tricky. While we conclude that the influence of vertical divergence is non-negligible, it has been ignored in some previous airborne flux studies (e.g., Vaughan et al.2016; Hannun et al.2020; Vaughan et al.2021; Drysdale et al.2022). While the flux is scattered in each vertical interval in our divergence calculation, we first bootstrap the flux observations and calculate the uncertainty in the correction factor (σC) to 40 %. As we see a significant difference in the vertical correction factor for racetrack measurements versus a selected subset of flux observations, we tentatively set the uncertainty of C to 100 % in order to account for the case of no vertical divergence. Besides, we account for a 30 % uncertainty in the PBL heights.

We propagate the total uncertainty from each component using Eq. (13), and the distribution of the total uncertainty of the 500 m average NOx flux is shown in Fig. 3b. The average uncertainty is 60 %, and the interquartile of the total uncertainty is 48 % and 68 %. The random error and the vertical divergence correction dominate the uncertainty, and the uncertainty is consistent with previous studies (Wolfe et al.2018; Vaughan et al.2016).

4 Component flux disaggregation

The overview of observed fluxes across seven flights over San Joaquin Valley is illustrated in Fig. 4. It shows a distinct spatial heterogeneity (Fig. 4a). For instance, high NOx flux signals are detected when the aircraft was flying above highway 99 between Bakersfield and Visalia. The transect of cities, such as Fresno, captures a substantial enhancement of NOx fluxes. Figure 4b exhibits the distribution of airborne fluxes. In total, 90 % of the fluxes are positive, demonstrating that our airborne flux measurements are capable of detecting NOx emissions over the study domain. We attribute the remaining 10 % of negative fluxes to the uncertainties in the flux calculation. The distribution of observed fluxes is right-skewed; the mean and median observed flux over the SJV is 0.37 and 0.25 mg N m−2 h−1, respectively. The interquartile range of the flux is 0.11 and 0.49 mg N m−2 h−1. In total, 1.2 % of extremely high fluxes exceeding 2 mg N m−2 h−1 represents the long tail in the flux distribution, which is, like the negative fluxes, most likely caused by the incomplete sampling of the spectrum of eddies driving the fluxes.

Figure 4(a) The map of observed airborne fluxes over seven flights over the San Joaquin Valley. If the segments overlap each other, then the average flux is calculated. The gray shading represents the coverage of 90 % footprint extents for all flux observations. Background map credits: Esri, HERE, Garmin, NGA, USGS, and NPS. (b) The distribution of full data sets of observed airborne NOx fluxes is shown.

As discussed in Sect. 3.3, we then calculate the footprint for each flux observation during the RECAP-CA field campaign. Figure 4a shows the 90 % footprint extent in gray. Figure S8 shows that the 90 % extent for the calculated footprints ranges from 0.16 to 12 km, with a mean extent of 2.8 km. The KL04-2D footprint algorithm has been applied to an airborne flux analysis over London, and in that study, the 90 % footprint extents range from 3 to 12 km from the measurement (Vaughan et al.2021). While the largest footprint extent is comparable with those from Vaughan et al. (2021), our calculated footprints mostly have a smaller extent, as 62 % of the footprint extents are within 3 km of the aircraft flight track. We attribute the small footprints to the stagnant weather conditions and weaker horizontal wind advection compared to London. The mean wind speed is 2.9 m s−1 for full observation data sets and 2.4 m s−1 for those data points with footprint extents less than 3 km. The largest footprint extent corresponds to observations at the foothills, due to the higher altitude relative to the boundary layer height and stronger horizontal wind advection.

The region covered by the footprints is composed of mixed land cover types. We use the 2018 United States Department of Agriculture (USDA) CropScape database (, last access: 24 August 2023) to describe the land cover types. The resolution has been degraded from the native 30 m resolution to 500 m. For each grid, the land cover type is assigned if a land type makes up more than 50 % of the 500 m grid cell. We generalize a “soil” land cover type if the land cover type is identified as being either croplands or grassland. The grids classified as “developed” in CropScape are dominated by anthropogenic activities, including transportation and fuel combustion. We overlay the national highway network and categorize the grids containing highways as “highway” land types. The remaining grids are classified as “urban”, and they correspond to the area with heavy populations in the absence of highways. The distinction between a “highway” and an “urban” land type is utilized to address on-road mobile sources. In total, 37 % of the flux observations include the highway land type in the 90 % footprint extent, 23 % of the observations include the urban land type, and 96 % of the observations include the cultivated soil land type.

To disentangle the flux emanating from different land cover types, we apply the Disaggregation combining Footprint analysis and Multivariate Regression (DFMR) methodology described in Hutjes et al. (2010). The observed fluxes are treated as the weighted sum of component fluxes from each land cover type as follows:

(14) F obs = k = 1 3 w k F k ,

where k1 to k3 denote the highway, urban, and soil land types, wk is the fractional area within the 90 % footprint contour, and Fk values are the corresponding component fluxes. The multi-linear regression is applied to observations from all flights, which consist of 4391 data points. We perform the Monte Carlo simulation to identify the uncertainty in the multi-linear regression due to the flux uncertainty. The resulting statistical uncertainty is shown in Fig. 5. The highway land type yields the highest flux of 0.96 mg N m−2 h−1, with a standard deviation of 0.04 mg N m−2 h−1. The areas classified as urban land type exhibit a flux of 0.43 (±0.02) mg N m−2 h−1, which is  50 % of the highway flux. Most likely, the fluxes from highway land types are even higher than 0.96 mg N m−2 h−1. Note that the land type map is at a 500 m spatial scale, so the grid classified as highway indeed includes both the highway and areas near the highway. If, for example, the highway is only 10 % of the true area of the land cover pixel, then the fluxes on the highway could be as much as 10 times larger. The cultivated soil land type flux of 0.30 (±0.01) mg N m−2 h−1 is large. It is about one-quarter the magnitude of the highway flux and half that of the urban flux. As the total area of soil pixels is much larger than the area of highway or urban pixels integrated across the SJV, cultivated soil NOx emissions are a major factor.

Figure 5Bootstrapped statistical results of multi-linear regression to resolve component fluxes from the highway, urban, and cultivated soil land types. Each bar represents the average component fluxes from each land type, and the black line shows the standard deviation.


5 Calculation of NOx emissions map using airborne NOx fluxes

While these separate component fluxes emphasize the distinction between individual land types at the spatial resolution of the land cover (500 m), we utilize the NOx fluxes to yield an estimate of NOx emissions at 4 km. For each 4 km grid, we collect the observed fluxes for which 90 % of the footprint overlaps with this grid area and define the weight rk as the fractional area that the footprint covers. The emissions (in units of mg N m−2 h−1) is calculated by the weighted average of the flux (Eq. 15). Only grids measured by at least five flux observations are considered in order to focus our attention on those pixels for which we have a statistically representative sample of the emissions.

(15) Emis i = k = 1 n 5 r k F k k = 1 n 5 r k

The emissions are calculated based on the observations from six flights during weekdays (Fig. 6a). The largest reported weekday emissions were on 3 June, when the median emissions were 0.39 mg N m−2 h−1. The lowest weekday emissions were observed on 15 June, with the median emissions of 0.14 mg N m−2 h−1. The large daily variation observed in estimated emissions during weekdays is partially due to the variation in flight routes and footprint coverage. This is illustrated by the daily estimated emissions map shown in Fig. S10.

As the emissions inventories make a distinction between weekdays and weekends and do not account for the daily variation on different weekdays, we average over the six weekday flights to yield the best estimate of emissions maps over the San Joaquin Valley derived from flux measurements (Fig. 6b). The median estimated weekday NOx emissions over the study domain is 0.26 mg N m−2 h−1, with an interquartile range of 0.14 and 0.46 mg N m−2 h−1. The observed emissions map describes high NOx emissions in the cities of Bakersfield (35.3 N, 119 W) and Fresno (36.75 N, 119.8 W) and along highway 99.

Figure 6(a) The box-and-whisker plot of observed emissions for each flight, aligned according to the order of the flight days. The box represents the interquartile ranges of observed emissions, and the line represents the median emissions. The whiskers show the maximum and minimum values. (b) The spatial distribution of emissions at 4 km over SJV as derived from observed fluxes during weekdays. The patch color shows the observed NOx emissions. The edge color denotes the land cover type, with the grid cells covering highways in white and those covering urban regions in black, and the rest of the grid cells are categorized with cultivated soil land cover types in green. © OpenStreetMap contributors 2022. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.

5.1 Evaluation of anthropogenic NOx emissions inventories

First, we compare the observations to the inventory developed by the California Air Resources Board (CARB). The anthropogenic emissions of NOx consist of mobile sources, stationary sources, and other NOx emissions from miscellaneous processes such as residential fuel combustion and managed disposal. In the CARB inventory, the mobile sources are estimated from the EMission FACtors (EMFAC) model, v1.0.2 (CARB2021a), and the Off-Road mobile source emissions model (CARB2021b). The stationary sources are estimated based on the reported survey of facilities within each local jurisdiction and the emissions factors from the California Air Toxics Emission Factor (CATEF) database (CARB2021c). Hereinafter, we utilize “EMFAC” to represent anthropogenic vehicle-related NOx emissions used in the CARB inventory. An alternative anthropogenic emissions inventory is the Fuel-Based Inventory for Vehicle Emissions (FIVE), developed by McDonald et al. (2012) and updated by Harkins et al. (2021). Both emissions inventories are at 4 km spatial resolution.

To disentangle the contribution of different NOx emissions sources, we attribute emissions at grid cells covering either highway or urban regions to anthropogenic emissions from transportation and fuel combustion, and those at remaining grid cells are categorized as soil NOx emissions. For each grid cell categorized as being dominant in anthropogenic emissions, we then match the emissions inventories representing the weekday scenario to the same hour and grids of emissions derived from measured fluxes. The corresponding hour of the estimated emissions is rounded to the closest hour of the observation times. Figures 7 and S12 show the comparison of observed anthropogenic emissions against EMFAC and FIVE emissions inventories. Over urban regions, the mean and median observed RECAP-CA NOx emissions are 0.37 mg N m−2 h−1, and the interquartile range is 0.14 and 0.58 mg N m−2 h−1. Both EMFAC and FIVE yield a good agreement with our measurements; the mean urban NOx emissions are 0.40 and 0.43 mg N m−2 h−1. However, the median urban NOx emissions in these inventories is 24 % and 22 % lower than the observation, respectively. The estimated NOx emissions on grid cells covering highways are more scattered. The median estimated NOx emissions is 0.24 mg N m−2 h−1. It is lower than on urban grid cells due to spatial averaging and due to the fact that most of the highway length is outside of the urban regions. The distribution of observed RECAP-CA NOx emissions from the highway is right-skewed and characterized by an interquartile range of 0.14 and 0.47 mg N m−2 h−1. We also note that over highway 99, the RECAP-CA NOx emissions are a factor of 3 higher than average on grid cells near congestion, reflecting the variation in the emissions caused by real-time traffic conditions. Both EMFAC and FIVE provide lower NOx emissions over highway grids; the median NOx emissions are 37 % and 50 % of those from the RECAP-CA observations. The highway pixels include land cover that is mostly non-highway and typically soil. If soil N emissions are substantially larger than in these inventories, then it is possible that the measurements and bottom-up inventories for highways are in better agreement than indicated by the figure.

Figure 7The box-and-whisker plot of observed RECAP-CA anthropogenic NOx emissions from transportation and fuel combustion, as well as those from EMFAC and FIVE emissions inventories, separated by highway and urban land cover types. The box is the interquartile range with the line of the median value. The maximum and minimum emissions are shown by whiskers, and the mean emissions are shown in red dots.


5.2 Evaluation of soil NOx scheme

Soil NOx emissions are determined by biogeochemical processes, including soil microbe‐mediated nitrification and denitrification. Process‐based biogeochemical models have been developed to mechanistically represent soil NOx emissions by simulating nitrogen interactions in ecological systems, such as DeNitrification‐-DeComposition (DNDC; Li et al.1992, 1994; Guo et al.2020) and DAYCENT (Del Grosso et al.2000; Rasool et al.2019). However, these process-level models are not yet widely applied to chemical transport models, and the default model configuration uses empirical soil NOx schemes. The Model of Emissions of Gases and Aerosols from Nature v3 (MEGAN; Guenther et al.2012) is the most commonly used scheme and is used to predict soil NOx emissions in the CARB emissions inventory. It is gridded at a 4 km spatial scale and has hourly time steps. The Biogenic Emission Inventory System (BEIS) is the default scheme to estimate volatile organic compounds from vegetation and NO from soil, as developed by the United States Environmental Protection Agency (EPA). We obtain the hourly BEIS v2.14 soil NOx emissions at 4 km during the study period from the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem), using the same model configuration described in Kim et al. (2022). While soil NOx varies nonlinearly with meteorological conditions, soil conditions, and agricultural activities, both MEGAN and BEIS simplify the nonlinearity to an activity factor (γ), as a function of ambient temperature, leaf area index, and leaf age. A recently developed soil NOx scheme, the Berkeley–Dalhousie–Iowa Soil NO Parameterization (BDISNP; Hudman et al.2012; Sha et al.2021), is an intermediate complexity model and includes more details than the other two in order to be more faithful to direct measurements made at soils and to describe their seasonal and hourly variations. The BDISNP model includes parameters representing the effects of soil moisture, temperature, and soil nitrogen, including fertilizer. Using the same WRF-Chem setup described in Sha et al. (2021), we also calculate the BDISNP soil NOx emissions during the study period at the spatial resolution of 2 km and regrid them to 4 km.

Figure 8a illustrates the range of soil N emissions derived from RECAP-CA observations, as compared to these three different soil NOx schemes. The analysis of the observations exhibits median cultivated soil NOx emissions of 0.26 mg N m−2 h−1; the interquartile range of the inferred emissions is 0.14 and 0.45 mg N m−2 h−1. MEGAN and BEIS both have emissions an order of magnitude lower, with median soil NOx emissions of 0.008 and 0.011 mg N m−2 h−1, respectively. The BDISNP soil NOx scheme shows median soil NOx emissions of 0.14 mg N m−2 h−1. Figure 8b exhibits a point-by-point comparison of the observed RECAP-CA and the BDISNP soil NOx emissions, showing that there is a correspondence between the two, but the model is 2.2 times lower than the observations. Figure S13a and d show the spatial distribution of soil NOx emissions from observation and BDISNP scheme. Both show higher soil NOx emissions between 35.75 and 36.25 N.

Figure 8(a) The box-and-whisker plot of observed soil NOx emissions and parameterized soil NOx emissions from MEGAN, BEIS, and BDISNP schemes. The mean soil NOx emissions are shown in black dots. (b) The scatterplot of soil NOx emissions calculated from BDISNP scheme and from flux measurements. The dashed black is the least square linear fit.


A distinct characteristic of soil NOx emissions is its temperature dependence. For instance, Oikawa et al. (2015) identified unusually high soil NOx emissions in a high-temperature agricultural region based on in situ observations. The temperature-driven increase in soil NOx emissions raises concerns in the future warmer climate, resulting in a larger contribution to O3 pollution (Romer et al.2018). Here we leverage our flux observations to probe this temperature dependence, as shown in Fig. 9. We collect observed NOx emissions for each flight and select the subset of NOx emissions on grids categorized as cultivated soil land type. We also collect corresponding mean soil temperature from WRF-Chem and match them to observed NOx emissions both in time and space. A range of soil temperature between 295 and 304 K is observed. We apply the linear regression fit between soil temperature and estimated soil NOx emissions and show a positive temperature dependence of soil NOx emissions with the slope of 0.02 (±0.008) mg N m−2 h−1 K−1. We then bin observed soil NOx emissions to three soil temperature categories, each of which has 2 K intervals. The median soil NOx emissions increase from 0.22 to 0.29 mg N m−2 h−1, with the median soil temperature increasing from 296 to 300 K. As the response to soil temperature is incorporated in the BDISNP scheme, we also bin the BDISNP parameterized soil NOx emissions into the same soil temperature categories. Both the RECAP-CA-measured flux and the BDISNP-modeled soil NOx emissions exhibit an approximately 33 % increase over the range of soil temperature shown.

Figure 9The dependence of soil NOx emissions on soil temperature from both the flux measurements (gray) and BDISNP scheme (orange). Both observed and BDISNP soil NOx emissions are binned, based on mean soil temperature from WRF-Chem. Three soil temperature bins are described with 4 K intervals. The box-and-whisker plot shows the distribution, and the black dot shows the mean within each bin. The line connects the median soil NOx emissions across three bins.


5.3 Discussion of soil NOx emissions

Soil NOx emissions in California have been studied in field experiments. Matson et al. (1997) measured soil NOx emissions from nine dominant crop types in SJV and reported mean fluxes of 0.01–0.09 mg N m−2 h−1. They also reported a large variation in the measured NOx flux among crops and among different fields of the same crop; the highest measured NOx flux is 0.17 mg N m−2 h−1, due to the fertilizer application and soil moisture characteristics. Horwath and Burger (2013) observed an average flux of 0.05–0.28 mg N m−2 h−1 at middays during the summertime from five crops in California, and the highest NOx flux is >4 mg N m−2 h−1 in systems receiving large N inputs, which results in high concentrations of ammonium. Oikawa et al. (2015) observed soil NOx emissions in a high-temperature fertilized agricultural region of the Imperial Valley, CA, ranging between −0.02 and 3.2 mg N m−2 h−1. They also conducted control experiments to investigate the soil NOx emissions response to fertilization and irrigation. The highest soil NOx flux was reported ∼10 d after the fertilizer at the soil volumetric water content of 30 % and a soil temperature of ∼313 K.

The mean soil NOx flux, 0.32 mg N m−2 h−1, derived in our flux measurements is higher than the mean fluxes reported in Matson et al. (1997) and Horwath and Burger (2013); however, the range of estimated soil NOx flux is within those in Horwath and Burger (2013) and Oikawa et al. (2015). Fertilizer is likely the primary contributor to the higher mean soil NOx flux in our study. The RECAP-CA field campaign was conducted in June, right after the month of peak fertilizer use in SJV (Guo et al.2020). As shown in Oikawa et al. (2015), the soil NOx flux can increase up to 5-fold within 20 d of fertilizer. The higher mean soil NOx flux is also contributed by higher soil temperature. In our study, the mean soil temperature is 299 K, with a range between 295 and 304 K, whereas the observations in Horwath and Burger (2013) and Oikawa et al. (2015) spread over a wider range of soil temperature (288–315 K). Consistent with our study, the temperature dependence of the soil NOx emissions is observed in these field experiments. Horwath and Burger (2013) reported a 2.5–3.5-fold increase in NOx fluxes with a 10 C increase in soil temperature. Oikawa et al. (2015) showed that the temperature dependence of soil NOx emissions is nonlinear; a steeper increase in soil NOx emissions was observed with the soil temperature exceeding 295 K.

It is worth noting the limitation of estimated soil NOx emissions in our study. First of all, we are unable to investigate the dependence of soil NOx emissions on meteorological drivers other than soil temperature, such as soil moisture, as modeled soil moisture presents a very small variation during the field campaign. Second, as our measurements only cover limited cropland areas in SJV over a short time period, and it is around the time of fertilizer use, we cannot scale the estimated soil NOx emissions to the whole year or to the total cropland areas in California. Therefore, we cannot directly compare our estimate of soil NOx emissions against other studies reporting soil NOx emissions on an annual basis or on a larger spatial scale. Last, in the absence of ozone and PM2.5 observations, we cannot investigate the impact of soil NOx emissions on air quality. However, as the SJV is in the NOx-limited regime (Pusede et al.2014), we expect that a model that captures the soil NOx more accurately will produce higher ozone. Future work is needed to further advance our understanding of soil NOx emissions and their role in urban and rural air pollution.

6 Conclusions

We performed airborne NOx flux measurements during RECAP-CA field campaign over the San Joaquin Valley. Seven flights were made over the SJV in June 2021. When combined with footprint and land cover information, we resolved the spatial heterogeneity in landscape flux. The component fluxes are estimated based on the multi-linear regression and exhibit statistically significant differences. The component fluxes are the highest from highways at 0.96 mg N m−2 h−1. Cultivated soil land types emit a non-negligible flux of 0.30 mg N m−2 h−1. The airborne flux observations are projected to a 4 km grid spacing to yield an estimated emissions map over the SJV. We utilize this map to evaluate emissions inventories commonly used in photochemical modeling. The anthropogenic emissions inventories, EMFAC and FIVE, agree well with the estimated mean NOx emissions over urban regions. However, the widely used, but not biogeochemical-process-based, models for soil NOx emissions underestimate emissions by an order of magnitude or more in the SJV, leading to a poor assessment of the relative roles of mobile and agriculture sources of NOx in the region. The BDISNP model, as adapted by Sha et al. (2021), results in a better comparison with the observations. Even though it is still lower by a factor of 2, we show that it yields a similar spatial pattern and soil temperature dependence, as observed. Variations in this model are embedded in CMAQ (Rasool et al.2019) and GEOS-CHEM (Wang et al.2021) and have been implemented in WRF-CHEM by Sha et al. (2021). Studies in which soil NOx is potentially important should make use of these codes, all of which are more consistent with observations at multiple scales.

Code and data availability

The measurement data from the RECAP-CA field champaign are available at (Chemical Sciences Laboratory2023), and the password is available upon request to the corresponding authors. The analysis codes for this study are available at (Zhu et al.2023).


The supplement related to this article is available online at:

Author contributions

RCC and AHG supervised the research. BP, EYP, BCS, PW CA, AB, JHS, RCC, and AHG participated in the field campaign. BP and PW conducted the NOx measurements. ST, HZ, and JW provided model-simulated BDISNP soil NOx emissions. QZ performed the analysis, with contributions from BP, EYP, and CMN. QZ prepared the paper. All authors have reviewed and edited the paper.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Atmospheric Chemistry and Physics. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors acknowledge the financial support of the California Air Resources Board and the Presidential Early Career Award for Scientists and Engineers (PECASE; from Brian McDonald) for the RECAP-CA field campaign. Qindan Zhu has been supported by the NOAA Climate and Global Change Postdoc Fellowship, and Eva Y. Pfannerstill has been supported by a Feodor Lynen Fellowship from the Alexander von Humboldt Foundation. Jun Wang and Huanxin Zhang acknowledge the support of NASA Atmospheric Composition and Modeling Program. We thank Dennis Baldocchi, Glenn Wolfe, and Erin Delaria for their help in calculating vertical divergence. We extend our gratitude to Brian McDonald, Rebecca Schwantes, and Siyuan Wang for engaging in discussions at project meetings. We appreciate use of the emissions inventories provided by Modeling and Meteorology Branch at CARB and NOAA Chemical Sciences Laboratory. We acknowledge the help from pilots, Bryce Kujat and George Loudakis, during the RECAP-CA field campaign.

Financial support

The RECAP-CA field campaign has been supported by the California Air Resources Board (grant nos. 20AQP012 and 20RD003) and PECASE. Jun Wang and Huanxin Zhang have been supported by ACMAP (grant no. 80NSSC19K0950).

Review statement

This paper was edited by Thomas Karl and reviewed by two anonymous referees.


Almaraz, M., Bai, E., Wang, C., Trousdell, J., Conley, S., Faloona, I., and Houlton, B. Z.: Agriculture is a major source of NOx pollution in California, Science Advances, 4, eaao3477,, 2018. a

American Lung Association: State of the Air Report, (last access: 22 June 2022), 2020. a

Andreae, M. O. and Schimel, D. S.: Exchange of trace gases between terrestrial ecosystems and the atmosphere, Plant Growth Regul., 10, 383–384,, 1990. a

CARB: Heavy-Duty Engine and Vehicle Omnibus Regulation and Associated Amendments, (last access: 22 June 2022), 2016. a

CARB: EMission FACtor (EMFAC), (last access: 22 June 2022), 2021a. a

CARB: OFFROAD2021 model, work/programs/mobile-source-emissions-inventory/msei-road-documentation-0 (last access: 22 June 2022), 2021b. a

CARB: Stationary Emissions by Air District, (last access: 22 June 2022), 2021c. a

Chemical Sciences Laboratory (CSL): RECAP-CA field champaign,, last access: 24 August 2023. a

Day, D., Wooldridge, P., Dillon, M., Thornton, J., and Cohen, R.: A thermal dissociation laser-induced fluorescence instrument for in situ detection of NO2, peroxy nitrates, alkyl nitrates, and HNO3, J. Geophys. Res.-Atmos., 107, ACH 4-1–ACH 4-14, 2002. a

De Gouw, J. A., Parrish, D. D., Frost, G. J., and Trainer, M.: Reduced emissions of CO2, NOx, and SO2 from US power plants owing to switch from coal to natural gas with combined cycle technology, Earth's Future, 2, 75–82, 2014. a

Del Grosso, S., Parton, W., Mosier, A., Ojima, D., Kulmala, A., and Phongpan, S.: General model for N2O and N2 gas emissions from soils due to dentrification, Global Biogeochem. Cy., 14, 1045–1060, 2000. a

Drysdale, W. S., Vaughan, A. R., Squires, F. A., Cliff, S. J., Metzger, S., Durden, D., Pingintha-Durden, N., Helfter, C., Nemitz, E., Grimmond, C. S. B., Barlow, J., Beevers, S., Stewart, G., Dajnak, D., Purvis, R. M., and Lee, J. D.: Eddy covariance measurements highlight sources of nitrogen oxide emissions missing from inventories for central London, Atmos. Chem. Phys., 22, 9413–9433,, 2022. a

EPA: Air Pollutant Emissions Trends Data, US EPA [data set], (last access: 28 October 2021), 2016. a

Gu, D., Guenther, A. B., Shilling, J. E., Yu, H., Huang, M., Zhao, C., Yang, Q., Martin, S. T., Artaxo, P., Kim, S., Seco, R., Stavrakou, T., Longo, K. M., Tóta, J., de Souza, R. A. F., Vega, O., Liu, Y., Shrivastava, M., Alves, E. G., Santos, F. C., Leng, G., and Hu, Z.: Airborne observations reveal elevational gradient in tropical forest isoprene emissions, Nat. Commun., 8, 1–7, 2017. a

Guenther, A. B., Jiang, X., Heald, C. L., Sakulyanontvittaya, T., Duhl, T., Emmons, L. K., and Wang, X.: The Model of Emissions of Gases and Aerosols from Nature version 2.1 (MEGAN2.1): an extended and updated framework for modeling biogenic emissions, Geosci. Model Dev., 5, 1471–1492,, 2012. a

Guo, L., Chen, J., Luo, D., Liu, S., Lee, H. J., Motallebi, N., Fong, A., Deng, J., Rasool, Q. Z., Avise, J. C., Kuwayama, T., Croes, B. E., and FitzGibbon, M.: Assessment of nitrogen oxide emissions and San Joaquin Valley PM2.5 impacts from soils in California, J. Geophys. Res.-Atmos., 126, e2020JD033304,, 2020. a, b, c

Hakeem, K. R., Sabir, M., Ozturk, M., Akhtar, M., and Ibrahim, F. H.: Nitrate and nitrogen oxides: sources, health effects and their remediation, Rev. Environ. Contam. T., 242, 183–217, 2016. a

Hannun, R. A., Wolfe, G. M., Kawa, S. R., Hanisco, T. F., Newman, P. A., Alfieri, J. G., Barrick, J., Clark, K. L., DiGangi, J. P., Diskin, G. S., King, J., Kustas, W. P., Mitra, B., Noormets, A., Nowak, J. B., Thornhill, K. L., and Vargas, R.: Spatial heterogeneity in CO2, CH4, and energy fluxes: Insights from airborne eddy covariance measurements over the Mid-Atlantic region, Environ. Res. Lett., 15, 035008,, 2020. a, b

Harkins, C., McDonald, B. C., Henze, D. K., and Wiedinmyer, C.: A fuel-based method for updating mobile source emissions during the COVID-19 pandemic, Environ. Res. Lett., 16, 065018,, 2021. a

Horwath, W. R. and Burger, M.: Assessment of NOx Emissions from Soil in California Cropping Systems, Tech. Rep. Contract No. 09–329, California Environmental Protection Agency, Air Resources Board, Sacramento, CA, (last access: 24 August 2023), 2013. a, b, c, d, e

Hudman, R. C., Moore, N. E., Mebust, A. K., Martin, R. V., Russell, A. R., Valin, L. C., and Cohen, R. C.: Steps towards a mechanistic model of global soil nitric oxide emissions: implementation and space based-constraints, Atmos. Chem. Phys., 12, 7779–7795,, 2012. a

Hutjes, R., Vellinga, O., Gioli, B., and Miglietta, F.: Dis-aggregation of airborne flux measurements using footprint analysis, Agr. Forest Meteorol., 150, 966–983, 2010. a

Kampa, M. and Castanas, E.: Human health effects of air pollution, Environ. Pollut., 151, 362–367, 2008. a

Karl, T., Apel, E., Hodzic, A., Riemer, D. D., Blake, D. R., and Wiedinmyer, C.: Emissions of volatile organic compounds inferred from airborne flux measurements over a megacity, Atmos. Chem. Phys., 9, 271–285,, 2009. a

Karl, T., Misztal, P., Jonsson, H., Shertz, S., Goldstein, A., and Guenther, A.: Airborne flux measurements of BVOCs above Californian oak forests: Experimental investigation of surface and entrainment fluxes, OH densities, and Damköhler numbers, J. Atmos. Sci., 70, 3277–3287, 2013. a, b, c

Kaser, L., Karl, T., Yuan, B., Mauldin III, R., Cantrell, C., Guenther, A. B., Patton, E., Weinheimer, A. J., Knote, C., Orlando, J., Emmons, L., Apel, E., Hornbrook, R., Shertz, S., Ullmann, K., Hall, S., Graus, M., de Gouw, J., Zhou, X., and Ye, C.: Chemistry-turbulence interactions and mesoscale variability influence the cleansing efficiency of the atmosphere, Geophys. Res. Lett., 42, 10894–10903, 2015. a

Kim, S.-W., McDonald, B. C., Seo, S., Kim, K.-M., and Trainer, M.: Understanding the paths of surface ozone abatement in the Los Angeles Basin, J. Geophys. Res.-Atmos., 127, e2021JD035606,, 2022. a

Kljun, N., Calanca, P., Rotach, M., and Schmid, H.: A simple parameterisation for flux footprint predictions, Bound.-Lay. Meteorol., 112, 503–523, 2004. a, b

Langford, B., Acton, W., Ammann, C., Valach, A., and Nemitz, E.: Eddy-covariance data with low signal-to-noise ratio: time-lag determination, uncertainties and limit of detection, Atmos. Meas. Tech., 8, 4197–4213,, 2015. a

Lelieveld, J., Evans, J. S., Fnais, M., Giannadaki, D., and Pozzer, A.: The contribution of outdoor air pollution sources to premature mortality on a global scale, Nature, 525, 367–371, 2015. a

Lenschow, D., Mann, J., and Kristensen, L.: How long is long enough when measuring fluxes and other turbulence statistics?, J. Atmos. Ocean. Tech., 11, 661–673, 1994. a, b

Li, C., Frolking, S., and Frolking, T. A.: A model of nitrous oxide evolution from soil driven by rainfall events: 1. Model structure and sensitivity, J. Geophys. Res.-Atmos., 97, 9759–9776, 1992. a

Li, C., Frolking, S., and Harriss, R.: Modeling carbon biogeochemistry in agricultural soils, Global Biogeochem. Cy., 8, 237–254, 1994. a

Matson, P., Firestone, M., Herman, D., Billow, T., Kiang, N., Benning, T., and Burns, J.: Agricultural Systems in the San Joaquin Valley: Development of Emissions Estimates for Nitrogen Oxides, Tech. Rep. Contract No. 94–732, California Environmental Protection Agency, Air Resources Board, Sacramento, CA, (last access: 24 August 2023), 1997. a, b

Mauder, M., Desjardins, R. L., and MacPherson, I.: Scale analysis of airborne flux measurements over heterogeneous terrain in a boreal ecosystem, J. Geophys. Res., 112, D13112,, 2007. a, b

McDonald, B. C., Dallmann, T. R., Martin, E. W., and Harley, R. A.: Long-term trends in nitrogen oxide emissions from motor vehicles at national, state, and air basin scales, J. Geophys. Res., 117, D00V18,, 2012. a

Metzger, S., Junkermann, W., Mauder, M., Beyrich, F., Butterbach-Bahl, K., Schmid, H. P., and Foken, T.: Eddy-covariance flux measurements with a weight-shift microlight aircraft, Atmos. Meas. Tech., 5, 1699–1717,, 2012. a

Misztal, P. K., Karl, T., Weber, R., Jonsson, H. H., Guenther, A. B., and Goldstein, A. H.: Airborne flux measurements of biogenic isoprene over California, Atmos. Chem. Phys., 14, 10631–10647,, 2014. a

Nussbaumer, C. M., Place, B. K., Zhu, Q., Pfannerstill, E. Y., Wooldridge, P., Schulze, B. C., Arata, C., Ward, R., Bucholtz, A., Seinfeld, J. H., Goldstein, A. H., and Cohen, R. C.: Measurement report: Airborne measurements of NOx fluxes over Los Angeles during the RECAP-CA 2021 campaign, EGUsphere [preprint],, 2023. a

Oikawa, P., Ge, C., Wang, J., Eberwein, J., Liang, L., Allsman, L., Grantz, D., and Jenerette, G.: Unusually high soil nitrogen oxide emissions influence air quality in a high-temperature agricultural region, Nat. Commun., 6, 1–10, 2015. a, b, c, d, e, f

Pfannerstill, E. Y., Arata, C., Zhu, Q., Schulze, B. C., Woods, R., Seinfeld, J. H., Bucholtz, A., Cohen, R. C., and Goldstein, A. H.: Volatile organic compound fluxes in the San Joaquin Valley – spatial distribution, source attribution, and inventory comparison, EGUsphere [preprint],, 2023. a, b, c

Phoenix, G. K., Hicks, W. K., Cinderby, S., Kuylenstierna, J. C., Stock, W. D., Dentener, F. J., Giller, K. E., Austin, A. T., Lefroy, R. D., Gimeno, B. S., Ashmore, M. R., and Ineson, P.: Atmospheric nitrogen deposition in world biodiversity hotspots: the need for a greater global perspective in assessing N deposition impacts, Glob. Change Biol., 12, 470–476, 2006. a

Pusede, S. E., Gentner, D. R., Wooldridge, P. J., Browne, E. C., Rollins, A. W., Min, K.-E., Russell, A. R., Thomas, J., Zhang, L., Brune, W. H., Henry, S. B., DiGangi, J. P., Keutsch, F. N., Harrold, S. A., Thornton, J. A., Beaver, M. R., St. Clair, J. M., Wennberg, P. O., Sanders, J., Ren, X., VandenBoer, T. C., Markovic, M. Z., Guha, A., Weber, R., Goldstein, A. H., and Cohen, R. C.: On the temperature dependence of organic reactivity, nitrogen oxides, ozone production, and the impact of emission controls in San Joaquin Valley, California, Atmos. Chem. Phys., 14, 3373–3395,, 2014. a

Rasool, Q. Z., Bash, J. O., and Cohan, D. S.: Mechanistic representation of soil nitrogen emissions in the Community Multiscale Air Quality (CMAQ) model v 5.1, Geosci. Model Dev., 12, 849–878,, 2019. a, b

Romer, P. S., Duffey, K. C., Wooldridge, P. J., Edgerton, E., Baumann, K., Feiner, P. A., Miller, D. O., Brune, W. H., Koss, A. R., de Gouw, J. A., Misztal, P. K., Goldstein, A. H., and Cohen, R. C.: Effects of temperature-dependent NOx emissions on continental ozone production, Atmos. Chem. Phys., 18, 2601–2614,, 2018. a

Sayres, D. S., Dobosy, R., Healy, C., Dumas, E., Kochendorfer, J., Munster, J., Wilkerson, J., Baker, B., and Anderson, J. G.: Arctic regional methane fluxes by ecotope as derived using eddy covariance from a low-flying aircraft, Atmos. Chem. Phys., 17, 8619–8633,, 2017. a

Schaller, C., Göckede, M., and Foken, T.: Flux calculation of short turbulent events – comparison of three methods, Atmos. Meas. Tech., 10, 869–880,, 2017. a

Schulze, B. C., Ward, R., Pfannerstill, E. Y., Zhu, Q., Arata, C., Place, B., Nussbaumer, C. M., Wooldridge, P., Woods, R., Bucholz, A., Cohen, R. C., Goldstein, A. H., Wennberg, P. O., and Seinfeld, J. H.: Methane emissions from dairy operations in California's San Joaquin Valley evaluated using airborne flux measurements, Environ. Sci. Technol., submitted, 2023. a

Sha, T., Ma, X., Zhang, H., Janechek, N., Wang, Y., Wang, Y., Castro García, L., Jenerette, G. D., and Wang, J.: Impacts of Soil NOx Emission on O3 Air Quality in Rural California, Environ. Sci. Technol., 55, 7113–7122, 2021. a, b, c, d, e

Sühring, M., Metzger, S., Xu, K., Durden, D., and Desai, A.: Trade-offs in flux disaggregation: a large-eddy simulation study, Bound.-Lay. Meteorol., 170, 69–93, 2019. a

Thornton, J. A., Wooldridge, P. J., and Cohen, R. C.: Atmospheric NO2: In situ laser-induced fluorescence detection at parts per trillion mixing ratios, Anal. Chem., 528–539, 2000. a

Torrence, C. and Compo, G. P.: A practical guide to wavelet analysis, B. Am. Meteorol. Soc., 79, 61–78, 1998. a

Tsai, J.-H., Huang, P.-H., and Chiang, H.-L.: Characteristics of volatile organic compounds from motorcycle exhaust emission during real-world driving, Atmos. Environ., 99, 215–226, 2014. a

Vaughan, A. R., Lee, J. D., Misztal, P. K., Metzger, S., Shaw, M. D., Lewis, A. C., Purvis, R. M., Carslaw, D. C., Goldstein, A. H., Hewitt, C. N., Davison, B., Beevers, S. D., and Karl, T. G.: Spatially resolved flux measurements of NOx from London suggest significantly higher emissions than predicted by inventories, Faraday Discuss., 189, 455–472, 2016. a, b

Vaughan, A. R., Lee, J. D., Metzger, S., Durden, D., Lewis, A. C., Shaw, M. D., Drysdale, W. S., Purvis, R. M., Davison, B., and Hewitt, C. N.: Spatially and temporally resolved measurements of NOx fluxes by airborne eddy covariance over Greater London, Atmos. Chem. Phys., 21, 15283–15298,, 2021. a, b, c, d

Wang, Y., Ge, C., Garcia, L. C., Jenerette, G. D., Oikawa, P. Y., and Wang, J.: Improved modelling of soil NOx emissions in a high temperature agricultural region: role of background emissions on NO2 trend over the US, Environ. Res. Lett., 16, 084061,, 2021. a

Wolfe, G., Hanisco, T., Arkinson, H., Bui, T., Crounse, J., Dean-Day, J., Goldstein, A., Guenther, A., Hall, S., Huey, G., Jacob, D. J., Karl, T., Kim, P. S., Liu, X., Marvin, M. R., Mikoviny, T., Misztal, P. K., Nguyen, T. B., Peischl, J., Pollack, I., Ryerson, T., St. Clair, J. M., Teng, A., Travis, K. R., Ullmann, K., Wennberg, P. O., and Wisthaler, A.: Quantifying sources and sinks of reactive gases in the lower atmosphere using airborne flux observations, Geophys. Res. Lett., 42, 8231–8240, 2015. a

Wolfe, G. M., Kawa, S. R., Hanisco, T. F., Hannun, R. A., Newman, P. A., Swanson, A., Bailey, S., Barrick, J., Thornhill, K. L., Diskin, G., DiGangi, J., Nowak, J. B., Sorenson, C., Bland, G., Yungel, J. K., and Swenson, C. A.: The NASA Carbon Airborne Flux Experiment (CARAFE): instrumentation and methodology, Atmos. Meas. Tech., 11, 1757–1776,, 2018. a, b, c

Wooldridge, P. J., Perring, A. E., Bertram, T. H., Flocke, F. M., Roberts, J. M., Singh, H. B., Huey, L. G., Thornton, J. A., Wolfe, G. M., Murphy, J. G., Fry, J. L., Rollins, A. W., LaFranchi, B. W., and Cohen, R. C.: Total Peroxy Nitrates (ΣPNs) in the atmosphere: the Thermal Dissociation-Laser Induced Fluorescence (TD-LIF) technique and comparisons to speciated PAN measurements, Atmos. Meas. Tech., 3, 593–607,, 2010.  a

Yu, H., Guenther, A., Gu, D., Warneke, C., Geron, C., Goldstein, A., Graus, M., Karl, T., Kaser, L., Misztal, P., and Yuan, B.: Airborne measurements of isoprene and monoterpene emissions from southeastern US forests, Sci. Total Environ., 595, 149–158, 2017. a

Yuan, B., Kaser, L., Karl, T., Graus, M., Peischl, J., Campos, T. L., Shertz, S., Apel, E. C., Hornbrook, R. S., Hills, A., Gilman, J. B., Lerner, B. M., Warneke, C., Flocke, F. M., Ryerson, T. B., Guenther, A. B., and de Gouw, J. A.: Airborne flux measurements of methane and volatile organic compounds over the Haynesville and Marcellus shale gas production regions, J. Geophys. Res.-Atmos., 120, 6271–6289, 2015. a, b

Zhu, Q., Place, B., Pfannerstill, E. Y., Tong, S., Zhang, H., Wang, J., Nussbaumer, C. M., Wooldridge, P., Schulze, B. C., Arata, C., Bucholtz, A., Seinfeld, J. H., Goldstein, A. H., and Cohen, R. C.: qdzhu/FLUX: v1.0 qdzhu/FLUX: First release, Version 1.0, Zenodo [code],, 2023. a

Short summary
Nitrogen oxide (NOx) is a hazardous air pollutant, and it is the precursor of short-lived climate forcers like tropospheric ozone and aerosol particles. While NOx emissions from transportation has been strictly regulated, soil NOx emissions are overlooked. We use the airborne flux measurements to observe NOx emissions from highways and urban and cultivated soil land cover types. We show non-negligible soil NOx emissions, which are significantly underestimated in current model simulations.
Final-revised paper