Atmos. Chem. Phys., 23, 9669–9683, 2023
https://doi.org/10.5194/acp-23-9669-2023
© Author(s) 2023. This work is distributed under
the Creative Commons Attribution 4.0 License.
Research article
Direct observations of NO
x
emissions over the San
Joaquin Valley using airborne flux measurements
during RECAP-CA 2021 field campaign
Qindan Zhu
1,*
, Bryan Place
2,**
, Eva Y. Pfannerstill
3
, Sha Tong
4,5
, Huanxin Zhang
4
, Jun Wang
4
,
Clara M. Nussbaumer
1,6
, Paul Wooldridge
2
, Benjamin C. Schulze
7
, Caleb Arata
8
, Anthony Bucholtz
9
,
John H. Seinfeld
7
, Allen H. Goldstein
3,8
, and Ronald C. Cohen
1,2
1
Department of Earth and Planetary Sciences, University of California,
Berkeley, Berkeley, CA 94720, USA
2
Department of Chemistry, University of California, Berkeley, Berkeley, CA 94720, USA
3
Department of Environmental Science, Policy, and Management, University of California,
Berkeley, Berkeley, CA 94720, USA
4
Department of Chemical and Biochemical Engineering, Center for Global and Regional Environmental
Research, and Iowa Technology Institute, University of Iowa, Iowa City, IA 52242, USA
5
Collaborative Innovation Center on Forecast and Evaluation of Meteorological Disasters (CIC-FEMD), Key
Laboratory for Aerosol-Cloud-Precipitation of China Meteorological Administration, Nanjing University of
Information Science and Technology, Nanjing 210044, People’s Republic of China
6
Department of Atmospheric Chemistry, Max Planck Institute for Chemistry, Mainz 55128, Germany
7
Department of Environmental Science and Engineering, California Institute of Technology,
Pasadena, CA 91125, USA
8
Department of Civil and Environmental Engineering, University of California,
Berkeley, Berkeley, CA 94720, USA
9
Department of Meteorology, Naval Postgraduate School, Monterey, CA 93943, USA
*
now at: Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology,
Cambridge, MA, United States
**
now at: Office of Research and Development, U.S. Environmental Protection Agency,
Research Triangle Park, NC 27711, USA
Correspondence:
Qindan Zhu (qindan_zhu@berkeley.edu) and Ronald C. Cohen (rccohen@berkeley.edu)
Received: 3 January 2023 – Discussion started: 24 February 2023
Revised: 28 June 2023 – Accepted: 12 July 2023 – Published: 31 August 2023
Abstract.
Nitrogen oxides (NO
x
) are principle components of air pollution and serve as important ozone pre-
cursors. As the San Joaquin Valley (SJV) experiences some of the worst air quality in the United States, reducing
NO
x
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 character-
ized 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 NO
x
emissions
using flux observations are utilized to evaluate anthropogenic emissions inventories and soil NO
x
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 mea-
sured fluxes over urban regions. Three soil NO
x
schemes, including the MEGAN v3 (Model of Emissions of
Published by Copernicus Publications on behalf of the European Geosciences Union.
9670
Q. Zhu et al.: NO
x
flux
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 NO
x
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 NO
x
emissions
and BDISNP present a similar spatial pattern and temperature dependence. We conclude that soil NO
x
is a key
feature of the NO
x
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 (NO
x
≡
NO
+
NO
2
) are important trace
gases that affect both the gas and aerosol phases of tro-
pospheric chemistry. NO
x
regulates the concentrations of
the primary atmospheric oxidant, hydroxyl radicals (OH),
and serves as the catalyst for the formation of ozone (O
3
).
NO
x
also affects the formation of inorganic nitrate aerosol
through the production of nitric acid (HNO
3
) and organic ni-
trates (RONO
2
) and plays a role in secondary organic aerosol
(SOA) production. NO
x
, O
3
, and aerosol are all detrimental
to human health, triggering respiratory diseases (Kampa and
Castanas, 2008; Hakeem et al., 2016) and leading to prema-
ture death (Lelieveld et al., 2015).
NO
x
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 (EPA,
2016). Strict regulations have been implemented to con-
trol NO
x
emissions. Three-way catalysts have effectively re-
duced emissions from gasoline-powered passenger vehicles.
The application of emissions control systems on coal power
plants has reduced NO
x
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 NO
x
emissions by 2031 (CARB,
2016). The regulation of mobile sources leads to an increas-
ing importance of natural NO
x
sources, such as lightning and
soil emissions. Soil NO
x
is released as a byproduct of micro-
bial nitrification and denitrification (Andreae and Schimel,
1990). While the biogeochemistry of soil NO
x
emissions is
well established, this biogenic source involves a complex in-
teraction 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 NO
x
emissions (Phoenix et al., 2006).
Currently, the San Joaquin Valley (SJV) in California ex-
periences 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 par-
ticulate matter (American Lung Association, 2020). In order
to implement appropriate emissions control efforts, identify-
ing the contribution of different NO
x
emissions is particu-
larly important for the SJV, as it features a complex mixture
of emissions from fuel combustion and soil emissions asso-
ciated with agriculture. The contribution of soil NO
x
emis-
sions remains highly uncertain. While Guo et al. (2020) at-
tribute approximately 1.1 % of anthropogenic NO
x
emissions
in California to soil NO
x
, Almaraz et al. (2018) argued that
due to growing N fertilizer use, the SJV has soil NO
x
emis-
sions of 24 kg N ha
−
1
yr
−
1
, contributing to 20 %–51 % of the
NO
x
budget of the entire state of California. Similarly, Sha et
al. (2021) estimated that 40.1 % of the total NO
x
emissions
over California in July 2018 are from soils.
Airborne eddy covariance (EC) flux measurements pro-
vide a powerful tool to investigate the emissions strength of
atmospheric constituents at landscape scales. It has been ap-
plied to assess the surface exchanges of greenhouse house
gases (GHGs), including CO
2
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 NO
x
over
a megacity (Karl et al., 2009; Vaughan et al., 2021), vege-
tation (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). Com-
pared to the traditional EC measurements from instruments
mounted at a fixed location on a tower, wavelet-based air-
borne 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 NO
x
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 sep-
arately. We utilize continuous wavelet transformation to cal-
culate the NO
x
flux (Sect. 3). In conjunction with the foot-
print calculations and land classification, we explore the spa-
tial heterogeneity of NO
x
emissions and identify component
fluxes from the highway, urban, and soil land types (Sect. 4).
We also utilize the NO
x
emissions derived from flux mea-
surements to evaluate anthropogenic emissions inventories
and soil NO
x
schemes (Sect. 5).
Atmos. Chem. Phys., 23, 9669–9683, 2023
https://doi.org/10.5194/acp-23-9669-2023
Q. Zhu et al.: NO
x
flux
9671
2 Measurements
The airborne EC flux measurements were conducted on a
Twin Otter research aircraft operated by the Naval Postgrad-
uate 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 tem-
perature, barometric and dynamic pressures, wind direction
and wind speed, total airspeed, slip and attack angles, GPS
latitude, GPS longitude, GPS altitude, pitch, roll, and head-
ing. These measurements are at 10 Hz temporal resolution.
VOCs were measured at 10 Hz time resolution by Vocus pro-
ton transfer reaction time-of-flight mass spectrometer (Vo-
cus PTR-ToF-MS), as described in Pfannerstill et al. (2023).
Mixing ratios of NO
x
were measured at 5 Hz frequency, us-
ing a custom-built three-channel thermal dissociation–laser-
induced fluorescence (TD–LIF) instrument. The multipass
LIF cells, fluorescence collection, long-pass wavelength fil-
tering (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 imple-
mentation 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 NO
2
by laser-induced flu-
orescence, 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 NO
2
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 cover-
age 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 NO
2
flu-
orescence lifetime for time-gated photon counting in order to
reject prompt laser scatter. Instrument zeros were run using
ambient air scrubbed of NO
x
every 20 min in flight to cor-
rect for any background drift during the flights. In addition,
calibrations were performed in-flight every 60 min, using a
NO
2
in N
2
calibration cylinder (Praxair; 5.5 ppm; Certified
Standard grade) diluted with scrubbed air.
NO
2
was measured directly in the first channel, with the
sample passing only through a particle filter and a flow-
limiting orifice before the cell. NO
x
was measured in the
second by adding O
3
(generated with 184.5 nm light and a
flow of scrubbed and dried air) to convert NO to NO
2
before
detection. A 122 cm length of 0.4 cm internal diameter tub-
ing served as the O
3
+
NO reactor, providing 4 s of reaction
time before the orifice. The third channel was used to mea-
sure the sum of higher nitrogen oxides (e.g., organic nitrates
and nitric acid) by thermal dissociation to NO
2
with an inline
oven (
∼
500
◦
C) 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 NO
x
mea-
surements. The full observation data set breaks into segments
with continuous wind and NO
x
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 activ-
ity. 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 concentra-
tion, 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 mea-
surements 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 tran-
sit time of air through the TD–LIF instrument, we assume
that the lag time for each flight is constant. We use the me-
dian lag time from each flight for all segments collected on
the same day.
3.2 Continuous wavelet transformation
The continuous wavelet transformation (CWT) parameteri-
zation decomposes the time series (
x
(
t
)) into a range of fre-
quencies and represents it as the convolution of the time se-
https://doi.org/10.5194/acp-23-9669-2023
Atmos. Chem. Phys., 23, 9669–9683, 2023
9672
Q. Zhu et al.: NO
x
flux
ries with a wavelet function (Torrence and Compo, 1998).
W
(
a,b
)
=
−∞
∫
∞
x
(
t
)
ψ
∗
a,b
(
t
)d
t
(1)
ψ
∗
a,b
(
t
)
=
1
√
a
ψ
0
(
t
−
b
a
)
,
(2)
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 loca-
tion of the “parent” wavelet, and the scale parameter defines
the frequency. We use the Morlet wavelet as the “parent”
wavelet.
ψ
0
=
π
−
1
/
4
e
6
iη
e
−
η
2
/
2
(3)
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 incre-
ment 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
a
j
with the increment
δj
of 0.25 s. The smallest frequency scale is the Nyquist fre-
quency, which is twice the time resolution (0.4 s).
b
n
=
Nδt
(4)
a
j
=
a
0
×
2
jδj
(5)
For two simultaneous time series of NO
x
(
W
c
(
a,b
)) and
vertical wind speed (
W
w
(
a,b
)), we obtain the wavelet cross-
spectrum, following Eq. (6). The Morlet-wavelet-specific re-
construction factor
C
δ
is 0.776. We then sum up over the full
frequency scales to yield a time series of the flux (Eq. 7).
E
c
,
w
(
j
)
=
δt
C
δ
1
N
N
−
1
∑
n
=
0
[
W
c
(
a,b
)
·
W
∗
w
(
a,b
)
]
(6)
F
(
t
)
=
c
′
w
′
=
δt
C
δ
δj
N
N
−
1
∑
n
=
0
J
∑
j
=
0
[
W
c
(
a,b
)
·
W
∗
w
(
a,b
)
]
a
(
j
)
(7)
Figure 1 exhibits an example of CWT flux calculation.
Figure 1a shows the normalized NO
x
and vertical wind speed
in a straight segment of
∼
50 km. The normalization is real-
ized 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 uncer-
tainties 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 spec-
trum is then integrated over all frequencies to the time series
of the NO
x
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.
3.3 Footprint calculation
The footprint describes the contribution of surface regions
to the observed airborne flux. We use the KL04-2D param-
eterization (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 disper-
sion model (Kljun et al., 2004). Metzger et al. (2012) im-
plemented a Gaussian cross-wind distribution function to re-
solve the dispersion perpendicular to the main wind direc-
tion. The input parameters include the height of the mea-
surements, standard deviation of the horizontal and verti-
cal 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 cal-
culate 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 geo-
graphic coordinate space.
3.4 Filter out NO
x
fluxes impacted by the off-road
vehicle emissions
It is worth noting that croplands include not only soil NO
x
emissions but also the off-road vehicle emissions. Erro-
neously attributing the NO
x
from off-road vehicle emis-
sions to soil NO
x
emissions leads to a high bias. While
trimethylbenzene was observed during RECAP-CA field
campaign, Pfannerstill et al. (2023) presented the trimethyl-
benzene fluxes using the same algorithm described in
Sect. 3.2. The trimethylbenzene fluxes are interpolated to
match the NO
x
fluxes in time and are utilized as an in-
dicator 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 foot-
prints 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 in-
cludes highway and urban areas. No correlation between
trimethylbenzene flux and NO
x
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
,
Atmos. Chem. Phys., 23, 9669–9683, 2023
https://doi.org/10.5194/acp-23-9669-2023
Q. Zhu et al.: NO
x
flux
9673
Figure 1.
(a)
The variance in the normalized NO
x
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.
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 trimethyl-
benzene flux between 0.005 and 0.04 mg m
−
2
h
−
1
and con-
clude 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 heterogene-
ity. 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 measure-
ments, 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 in-
cluded to avoid emissions heterogeneity. As shown in Fig. 2,
we calculate the ratio of measurement heights relative to the
PBL height (
z/z
i
), and 98 % of selected fluxes are located
within 70 % of the PBL height and are divided into seven
bins of
z/z
i
with uniform width. We then perform a linear fit
for the binned median fluxes versus
z/z
i
to calculate the ver-
tical correction factor
(
C
=
slope
intercept
)
. This correction factor
is used to linearly extrapolated the fluxes at the measurement
height (
F
z
) to fluxes at the surface (
F
0
) (Eq. 8). After the ver-
tical divergence correction, the surface fluxes are on average
26 % higher than the fluxes at the measurement heights.
F
0
=
F
x
1
+
C
z
z
i
(8)
3.6 Data qualify control and uncertainty analysis
The flux detection limit does not only depend on the signal-
to-noise ratio of the NO
x
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 NO
x
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 er-
ror affecting the flux (
σ
NO
x
,
noise
) is defined as the standard
deviation of this noise-derived flux, and LoD is defined as
2
×
σ
NO
x
,
noise
(95th confidence level). Among 142 segments,
Fig. 3a shows the distribution of flux LoD among 142 seg-
ments. 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 uncer-
tainty from a variety of sources. We describe the systematic
errors and random errors, following the work of Wolfe et al.
https://doi.org/10.5194/acp-23-9669-2023
Atmos. Chem. Phys., 23, 9669–9683, 2023
9674
Q. Zhu et al.: NO
x
flux
Figure 2.
Vertical profiles of measured fluxes above croplands dur-
ing RECAP-CA field campaign binned by the ratio of measurement
height and PBL height (
z/z
i
). 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 rela-
tive height.
Figure 3.
(a)
The distribution of the segment-based NO
x
flux de-
tection limit (LoD).
(b)
The distribution of total uncertainty in the
NO
x
flux at 500 m spatial average.
(2018). Systematic errors arise from the undersampling of
high-frequency and low-frequency ranges. The CWT algo-
rithm 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).
SE
≤
2
.
2
(
z
z
i
)
0
.
5
z
i
L
(9)
z
and
L
are the measurement heights and the length of seg-
ments, respectively.
z
i
values are the boundary layer heights
from HRRR. We calculate the low-frequency error ranges
from 1 % to 5 %.
Random errors arise from the noise in the instrument
(RE
noise
) and the noise in turbulence sampling (RE
turb
),
which are calculated using Eqs. (10) and (11) (Wolfe et al.,
2018; Lenschow et al., 1994).
RE
noise
=
√
σ
2
NO
x
,
noise
σ
2
w
N
(10)
RE
turb
F
≤
1
.
75
(
z
z
i
)
0
.
25
(
z
i
L
)
0
.
5
(11)
z
,
L
, and
z
i
values are the same as Eq. (9), and
σ
2
w
is the vari-
ance of the vertical wind speed. Note that RE
noise
assumes
that the noise in each time step is uncorrelated; therefore, we
ignore the moving average step in the uncertainty calcula-
tion, 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 com-
paring the calculated fluxes using segment-specific and con-
stant 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 (RE
lag
).
While we believe that this error is unphysical and that a sin-
gle lag time is more appropriate, we include it to be conser-
vative in our estimate of the uncertainties.
Estimating the uncertainty caused by the correction of ver-
tical divergence is tricky. While we conclude that the influ-
ence of vertical divergence is non-negligible, it has been ig-
nored in some previous airborne flux studies (e.g., Vaughan
et al., 2016; Hannun et al., 2020; Vaughan et al., 2021; Drys-
dale 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 correc-
tion factor (
σ
C
) to 40 %. As we see a significant difference in
the vertical correction factor for racetrack measurements ver-
sus 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 % un-
certainty 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 NO
x
flux is shown in Fig. 3b. The av-
erage uncertainty is 60 %, and the interquartile of the total
uncertainty is 48 % and 68 %. The random error and the ver-
tical divergence correction dominate the uncertainty, and the
uncertainty is consistent with previous studies (Wolfe et al.,
2018; Vaughan et al., 2016).
σ
F
z
=
√
SE
2
+
RE
2
noise
+
RE
2
turb
+
RE
2
lag
(12)
σ
F
0
=
√
√
√
√
√
√
σ
2
F
z
(
1
+
C
z
z
i
)
2
+
σ
2
C
(
z
z
i
)
2
F
z
(
1
+
C
z
z
i
)
2
2
+
σ
2
z
i
(
Cz
z
2
i
)
2
F
z
(
1
+
C
z
z
i
)
2
2
(13)
Atmos. Chem. Phys., 23, 9669–9683, 2023
https://doi.org/10.5194/acp-23-9669-2023
Q. Zhu et al.: NO
x
flux
9675
4 Component flux disaggregation
The overview of observed fluxes across seven flights over
San Joaquin Valley is illustrated in Fig. 4. It shows a dis-
tinct spatial heterogeneity (Fig. 4a). For instance, high NO
x
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 NO
x
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
NO
x
emissions over the study domain. We attribute the re-
maining 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 interquar-
tile 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.
As discussed in Sect. 3.3, we then calculate the footprint
for each flux observation during the RECAP-CA field cam-
paign. Figure 4a shows the 90 % footprint extent in gray. Fig-
ure S8 shows that the 90 % extent for the calculated foot-
prints 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 foot-
print extent is comparable with those from Vaughan et al.
(2021), our calculated footprints mostly have a smaller ex-
tent, 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 Depart-
ment of Agriculture (USDA) CropScape database (https:
//nassgeodata.gmu.edu/CropScape/, 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 re-
maining grids are classified as “urban”, and they correspond
to the area with heavy populations in the absence of high-
ways. The distinction between a “highway” and an “urban”
land type is utilized to address on-road mobile sources. In to-
tal, 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 in-
clude the cultivated soil land type.
To disentangle the flux emanating from different land
cover types, we apply the Disaggregation combining Foot-
print analysis and Multivariate Regression (DFMR) method-
ology 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:
F
obs
=
3
∑
k
=
1
w
k
F
k
,
(14)
where
k
1
to
k
3
denote the highway, urban, and soil land
types,
w
k
is the fractional area within the 90 % footprint con-
tour, and
F
k
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 re-
sulting statistical uncertainty is shown in Fig. 5. The high-
way 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 ar-
eas 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 high-
way. 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 NO
x
emissions are a major
factor.
5 Calculation of NO
x
emissions map using airborne
NO
x
fluxes
While these separate component fluxes emphasize the dis-
tinction between individual land types at the spatial resolu-
tion of the land cover (500 m), we utilize the NO
x
fluxes to
yield an estimate of NO
x
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
r
k
as the fractional area that the footprint covers. The emissions
https://doi.org/10.5194/acp-23-9669-2023
Atmos. Chem. Phys., 23, 9669–9683, 2023
9676
Q. Zhu et al.: NO
x
flux
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 NO
x
fluxes is shown.
Figure 5.
Bootstrapped statistical results of multi-linear regression
to resolve component fluxes from the highway, urban, and culti-
vated soil land types. Each bar represents the average component
fluxes from each land type, and the black line shows the standard
deviation.
(in units of mg N m
−
2
h
−
1
) is calculated by the weighted av-
erage 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.
Emis
i
=
∑
n
≥
5
k
=
1
r
k
F
k
∑
n
≥
5
k
=
1
r
k
(15)
The emissions are calculated based on the observations
from six flights during weekdays (Fig. 6a). The largest re-
ported weekday emissions were on 3 June, when the me-
dian emissions were 0.39 mg N m
−
2
h
−
1
. The lowest week-
day 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 par-
tially due to the variation in flight routes and footprint cover-
age. 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 NO
x
emissions
over the study domain is 0.26 mg N m
−
2
h
−
1
, with an in-
terquartile range of 0.14 and 0.46 mg N m
−
2
h
−
1
. The ob-
served emissions map describes high NO
x
emissions in the
cities of Bakersfield (35.3
◦
N, 119
◦
W) and Fresno (36.75
◦
N,
119.8
◦
W) and along highway 99.
5.1 Evaluation of anthropogenic NO
x
emissions
inventories
First, we compare the observations to the inventory de-
veloped by the California Air Resources Board (CARB).
The anthropogenic emissions of NO
x
consist of mobile
sources, stationary sources, and other NO
x
emissions from
miscellaneous processes such as residential fuel combus-
tion and managed disposal. In the CARB inventory, the mo-
bile sources are estimated from the EMission FACtors (EM-
FAC) model, v1.0.2 (CARB, 2021a), and the Off-Road mo-
bile source emissions model (CARB, 2021b). The stationary
sources are estimated based on the reported survey of facil-
ities within each local jurisdiction and the emissions factors
from the California Air Toxics Emission Factor (CATEF)
database (CARB, 2021c). Hereinafter, we utilize “EMFAC”
to represent anthropogenic vehicle-related NO
x
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
Atmos. Chem. Phys., 23, 9669–9683, 2023
https://doi.org/10.5194/acp-23-9669-2023