Airborne observation with a low-cost hyperspectral instrument: retrieval of NO2 vertical column densities (VCDs) and the satellite sub-grid variability over industrial point sources

Abstract

High-spatial-resolution NO${}_{\mathrm{2}}$ vertical column densities (VCDs) were retrieved from airborne observations using the low-cost hyperspectral imaging sensor (HIS) at three industrial areas (i.e., Chungnam, Jecheon, and Pohang) in South Korea, where point sources (i.e., power plant, petrochemical complex, steel yard, and cement kiln) with significant NO${}_{\mathrm{2}}$ emissions are located. An innovative and versatile approach for NO${}_{\mathrm{2}}$ VCD retrieval, hereafter referred to as the modified wavelength pair (MWP) method, was developed to overcome the excessively variable radiometric and spectral characteristics of the HIS attributed to the absence of temperature control during the flight. The newly developed MWP method was designed to be insensitive to broadband spectral features, including the spectral dependency of surface and aerosol reflectivity, and can be applied to observations with relatively low spectral resolutions. Moreover, the MWP method can be implemented without requiring precise radiometric calibration of the instrument (i.e., HIS) by utilizing clean-pixel data for non-uniformity corrections and is also less sensitive to the optical properties of the instrument and offers computational cost competitiveness. In the experimental flights using the HIS, NO${}_{\mathrm{2}}$ plumes emitted from steel yards were particularly conspicuous among the various NO${}_{\mathrm{2}}$ point sources, with peak NO${}_{\mathrm{2}}$ VCDs of 2.0 DU (Dobson unit) at Chungnam and 1.8 DU at Pohang. Typical NO${}_{\mathrm{2}}$ VCD uncertainties ranged between 0.025–0.075 DU over the land surface and 0.10–0.15 DU over the ocean surface, and the discrepancy can be attributable to the lower signal-to-noise ratio over the ocean and higher sensitivity of the MWP method to surface reflectance uncertainties under low-albedo conditions. The NO${}_{\mathrm{2}}$ VCDs retrieved from the HIS with the MWP method showed a good correlation with the collocated Tropospheric Monitoring Instrument (TROPOMI) data ($r=\mathrm{0.73}$, mean absolute error equals 0.106 DU). However, the temporal disparities between the HIS frames and the TROPOMI overpass, their spatial mismatch, and their different observation geometries could limit the correlation. The comparison of TROPOMI and HIS NO${}_{\mathrm{2}}$ VCDs further demonstrated that the satellite sub-grid variability could be intensified near the point sources, with more than a 3-fold increase in HIS NO${}_{\mathrm{2}}$ VCD variability (e.g., difference between 25th and 75th quantiles) over the TROPOMI footprints with NO${}_{\mathrm{2}}$ VCD values exceeding 0.8 DU compared to footprints with NO${}_{\mathrm{2}}$ VCD values below 0.6 DU.

1 Introduction

Reactive nitrogen oxides (${\mathrm{NO}}_x$), which is a sum of nitrogen dioxide (NO${}_{\mathrm{2}}$) and nitrogen monoxide (NO), are one of the most important trace gases in a polluted atmosphere. ${\mathrm{NO}}_x$ has been widely known to have direct health effects (Jo et al., 2021; Kampa and Castanas, 2008; Song et al., 2023) and plays an essential role in secondary ozone (O${}_{\mathrm{3}}$) and aerosol formation (Chan et al., 2010; Seinfeld and Pandis, 2016; Sillman, 1999), both of which are also the major atmospheric pollutants. ${\mathrm{NO}}_x$ can be emitted from natural sources such as wildfires or lightning (Jaffe and Wigder, 2012; Stark et al., 1996), whereas anthropogenic emissions, mainly from the combustion of fossil fuels, are the main cause of the air quality deterioration in populated urban and industrial regions (Kim et al., 2013; Lamsal et al., 2013). ${\mathrm{NO}}_x$ is highly reactive in the polluted atmosphere and exhibits large spatial inhomogeneity and a short atmospheric lifetime (Beirle et al., 2011; Liu et al., 2016; Valin et al., 2013). Therefore, surface in situ measurements of speciated ${\mathrm{NO}}_x$ must be complemented by remote sensing observations with wider spatial coverage. NO${}_{\mathrm{2}}$ has a strong and highly structured absorption cross-section at the visible (VIS) band, and the column densities can be retrieved from passive hyperspectral observations and with various spectral analysis techniques such as differential optical absorption spectroscopy (DOAS; Platt and Stutz, 2008). Considering the balanced state of NO and NO${}_{\mathrm{2}}$ from ${\mathrm{NO}}_x$ titration and photolysis during the daytime (Sillman, 1999), observations of NO${}_{\mathrm{2}}$ concentrations can yield the general concentration of ${\mathrm{NO}}_x$ in the atmosphere.

Since the first global NO${}_{\mathrm{2}}$ vertical column density (VCD) observations from the push-broom-type satellite hyperspectral imaging sensor, the Global Ozone Monitoring Experiment (GOME; Burrows et al., 1997), numerous hyperspectral imagers such as Ozone Monitoring Instrument (OMI; Levelt et al., 2006), GOME-2 (Munro et al., 2016), and Tropospheric Monitoring Instrument (TROPOMI; Veefkind et al., 2012) were sequentially launched with gradual enhancement in its spatial resolution. For instance, TROPOMI, the latest sensor in low Earth orbit, has a footprint of $\mathrm{3.5}\times \mathrm{5.5}$ km${}^{\mathrm{2}}$ at the nadir, while the GOME had a footprint of $\mathrm{40}\times \mathrm{320}$ km${}^{\mathrm{2}}$. However, even up-to-date environmental satellite sensors such as TROPOMI show inherent limitations in spatial resolution to fully capture the NO${}_{\mathrm{2}}$ spatial inhomogeneities in urban and industrial areas (Judd et al., 2019; Park et al., 2022; Tang et al., 2021; Verhoelst et al., 2021).

The airborne observation of NO${}_{\mathrm{2}}$ VCDs from hyperspectral imagers has demonstrated its effectiveness in elucidating the spatial variability in NO${}_{\mathrm{2}}$ near emission hot spots (Judd et al., 2019; Tack et al., 2021). However, the highly resolved airborne observations of NO${}_{\mathrm{2}}$ VCDs were sparse and limited to intensive field campaigns. State-of-the-art hyperspectral imagers dedicated to observing column densities of trace gases, such as the Geostationary Trace gas and Aerosol Sensor Optimization (GeoTASO; Leitch et al., 2014), GeoCAPE Airborne Simulator (GCAS; Kowalewski and Janz, 2014), Airborne imaging DOAS instrument for Measurements of Atmospheric Pollution (AirMAP; Schönhardt et al., 2015), and Airborne Compact Atmospheric Mapper (ACAM; Kowalewski and Janz, 2009; Lamsal et al., 2017), have been widely deployed in intensive field campaigns. Although these hyperspectral radiometers provide highly accurate measurements owing to their stable spectral and radiometric characteristics and precise spectral resolution, these instruments are expensive to develop and require sophisticated and meticulous calibration and maintenance to maintain their performance.

Meanwhile, efforts have been made to retrieve NO${}_{\mathrm{2}}$ VCDs from airborne observations using hyperspectral imagers that were not originally designed for trace-gas observations. Popp et al. (2012) and Tack et al. (2017, 2021) used the Airborne Prism Experiment (APEX; Itten et al., 2008), which is a hyperspectral imager initially developed to observe land use–land cover (LULC; Tack et al., 2019), to retrieve NO${}_{\mathrm{2}}$ VCDs in various airborne experiments. The APEX instrument has a variable and coarse spectral resolution in the VIS band with full width at half maximum (FWHM) ranging from 1.5 to 3.0 nm. Nevertheless, the NO${}_{\mathrm{2}}$ VCDs were successfully retrieved by applying conventional DOAS fitting and air mass factor (AMF) calculations, with a good agreement to collocated ground-based DOAS instruments ($r=\mathrm{0.84}$) and TROPOMI ($r=\mathrm{0.94}$; APEX data upscaled to the TROPOMI footprints).

To enhance our understanding and obtain explicit depictions of atmospheric pollutants, improving the feasibility of airborne observations is important. The value of airborne NO${}_{\mathrm{2}}$ VCD observations can be further highlighted in East Asia, where ${\mathrm{NO}}_x$ emissions are still substantial, and tropospheric O${}_{\mathrm{3}}$ levels are increasing (Colombi et al., 2023; Li et al., 2020). However, there are challenges associated with airborne hyperspectral observations, including the development of sophisticated instruments and the need for precise calibration and maintenance. Therefore, in this study, we developed a versatile NO${}_{\mathrm{2}}$ VCD retrieval algorithm suitable for a broader range of hyperspectral instruments. The newly developed algorithm has been implemented for spectra obtained from experimental flights with a low-cost hyperspectral imaging sensor, and the associated retrieval uncertainties were estimated. Furthermore, retrieved NO${}_{\mathrm{2}}$ VCDs were compared with the spatially collocated TROPOMI data, focusing on the satellite sub-grid variability observed over industrial point sources in South Korea.

2 The low-cost hyperspectral imaging sensor (HIS)

The compact, low-cost hyperspectral imaging sensor (hereafter referred to as HIS) designed to observe from the ultraviolet (UV) to VIS bands (250–500 nm) in hyperspectral resolution was used in this study for airborne observations. Table 1 shows the detailed specifications of the HIS, manufactured by Headwall Photonics, Inc. The HIS has a two-dimensional 16-bit charge-coupled device (CCD), which supports a relatively large dynamic range of digital photon counts from 0 to ${\mathrm{2}}^{\mathrm{16}}$, with 1392 (spatial columns) $\times \mathrm{1040}$ (spectral rows) arrays as a detector. The CCD conversion rate to digital count as a response to incident photons is known to be nonlinear when the CCD counts exceed approximately 80 % of the saturation level (Nehir et al., 2019; Pulli et al., 2017; Wang et al., 2016), and thus the CCD counts should remain below 50 000 counts for practical use. Furthermore, the spectral rows at the edges of the CCD are designed not to be illuminated, yielding 945 effective spectral rows. A diffraction grating is used with a concave mirror to disperse light as a spectrum, and a 25 $\mathrm{\mu }\mathrm{m}$ width slit on the light entrance is designed to yield a slit function with 1.4 nm FWHM. The HIS has a 0.265 nm spectral binning interval on average and a moderate level of oversampling (approximately 5 times) accordingly.

Table 1

Specifications of the HIS used in this study.

Specifications	Values	Remarks
Wavelength range	250–500 nm	945 effective spectral pixels
Spectral sampling interval	0.265 nm	Spectral binning ($\pm \mathrm{2}$ pixels) applied for the retrieval
Spectral resolution	1.4 nm	FWHM
Spatial pixels	1392	–
FOV	13${}^{\circ }$	Varying with the focal length of the lens
Detector	Interline CCD	16-bit
Slit width	25 $\mathrm{\mu }\mathrm{m}$	–
$f$ number	$f/\mathrm{2.0}$	–
Lens EFL	28.3 mm	The same objective lens was used throughout the study

Despite these instrument specifications, the HIS shows almost no sensitivity beyond the UV-A region ($\mathit{\lambda }<\mathrm{320}$ nm). On the other hand, the HIS exhibits sufficient sensitivity at the VIS range (350–500 nm), where the spectral absorption characteristics of NO${}_{\mathrm{2}}$ are well distinguished, making it suitable for NO${}_{\mathrm{2}}$ VCD retrievals. For instance, Park et al. (2019) retrieved NO${}_{\mathrm{2}}$ VCDs from ground-based observations of the diffuse radiances using the identical instrument (i.e., HIS) and presented its compatibility for NO${}_{\mathrm{2}}$ VCD retrievals by demonstrating a high correlation ($r=\mathrm{0.84}$) with collocated Pandora observations. However, there are several vulnerabilities in the HIS that hinder its utilization for airborne observations, primarily due to low-quality optical and radiometric characteristics resulting from its compact configuration and lack of regular or meticulous maintenance.

First and foremost, the HIS is not equipped with a temperature control unit, which is crucial to ensure ground calibration to preserve its validity and to keep optical and radiometric properties stable. Exposure to fluctuating ambient temperatures during the flight can cause thermal shifts in the instrument's optics, resulting in spectral shifts or changes in the slit function. Furthermore, the sensitivity of the CCD pixels can vary, even from pixel to pixel, with temperature changes, leading to additional pixel non-uniformity. The HIS inherently exhibits relatively unstable spectral and radiometric characteristics compared to the sophisticated state-of-the-art hyperspectral imagers like GeoTASO or GCAS, which is another drawback. Additionally, the narrow field of view (FOV) of the HIS can be a critical defect when used for airborne observations, resulting in a narrow swath. Theoretically, the FOV can be enlarged by using an objective lens with a shorter focal length. However, the objective lens with an effective focal length (EFL) of 28.3 mm used in this study is already compact with an excessively short focal length (i.e., low $f$ number) and is not feasible to enlarge the field of view under the current instrument configuration.

3

NO${}_{\mathrm{2}}$ VCD retrieval algorithm

Most modern passive remote sensing techniques aimed to observe atmospheric trace gases (e.g., NO${}_{\mathrm{2}}$, O${}_{\mathrm{3}}$, formaldehyde) in UV–VIS bands via hyperspectral sensors adopt the DOAS technique, which necessarily requires sufficiently high spectral resolution and well-defined, stable optical characteristics (e.g., slit function). More specifically, the DOAS fitting technique is applied to the observed spectra from airborne or space-borne hyperspectral images to retrieve slant column densities (SCDs), and by dividing by the AMF calculated from the forward radiative transfer model (RTM) and the chemical transfer model (CTM), it yields VCDs. To follow the convention and retrieve VCDs accurately, precise spectral and radiometric characterizations before and after each research flight, as well as stabilization in terms of temperature and vibration during the flight, are necessary. However, the demand for precise calibration and the need for ancillary components, such as temperature control units, pose obstacles that impede the feasibility of airborne hyperspectral observation. Therefore, we have developed an innovative and versatile approach for NO${}_{\mathrm{2}}$ VCD retrieval suitable for applying to spectra obtained from various hyperspectral instruments, hereafter referred to as the modified wavelength pair (MWP) method, inspired by the wavelength pair methods previously developed to retrieve columnar concentrations of atmospheric O${}_{\mathrm{3}}$ and NO${}_{\mathrm{2}}$ from ground-based instruments such as the Dobson spectrophotometer (Dobson, 1957a, b) and the Brewer spectrophotometer (Brewer and Kerr, 1973; Brewer et al., 1973). Figure 1 illustrates the overall flowchart outlining the NO${}_{\mathrm{2}}$ VCD retrieval from the airborne HIS spectra, including the implementation of the newly developed MWP method. Each step of the flowchart is elaborated in the following sections.

Figure 1

Flowchart of the NO${}_{\mathrm{2}}$ VCD retrieval from the airborne HIS spectra using the MWP method.

[Figure omitted. See PDF]

3.1 Theoretical basis of the modified wavelength pair (MWP) method

A total of six wavelength pairs have been selected within the fitting window ranging from 410–450 nm, where the NO${}_{\mathrm{2}}$ absorption exhibits significant high-frequency spectral variability within the VIS band. These pairs are composed of wavelengths that exhibit weak and strong NO${}_{\mathrm{2}}$ absorption and are classified into two sub-types – Type_A and Type_B – with three wavelength pairs allocated to each type (Table 2). In Type_A pairs, the shorter wavelength ${\mathit{\lambda }}_{\mathrm{1}}$ shows stronger sensitivity to NO${}_{\mathrm{2}}$ concentration, while the longer wavelength ${\mathit{\lambda }}_{\mathrm{2}}$ exhibits relatively weaker NO${}_{\mathrm{2}}$ sensitivity (Fig. 2). In contrast, Type_B pairs show the opposite pattern, with weaker NO${}_{\mathrm{2}}$ sensitivity at the shorter wavelength (${\mathit{\lambda }}_{\mathrm{1}}$) and stronger sensitivity at the longer wavelength (${\mathit{\lambda }}_{\mathrm{2}}$). The basic concept of using the wavelength pair is that the radiance ratio ($R=I_{{\mathit{\lambda }}_{\mathrm{1}}}/I_{{\mathit{\lambda }}_{\mathrm{2}}}$) calculated from the spectral radiance at two different wavelengths of a wavelength pair, ${\mathit{\lambda }}_{\mathrm{1}}$ and ${\mathit{\lambda }}_{\mathrm{2}}$ (where ${\mathit{\lambda }}_{\mathrm{1}}<{\mathit{\lambda }}_{\mathrm{2}}$), is expected to be pseudo-linearly dependent on the NO${}_{\mathrm{2}}$ VCD. Therefore, NO${}_{\mathrm{2}}$ VCDs can be retrieved by comparing the observed $R$ value with the $R$ values simulated from the RTM for a range of possible NO${}_{\mathrm{2}}$ VCDs (and controlling other possible artifacts).

Table 2

Wavelength pairs used in the NO${}_{\mathrm{2}}$ VCD retrieval algorithm (MWP method).

MWP	Wavelength	Wavelength 1	Wavelength 2
type	pair index	(${\mathit{\lambda }}_{\mathrm{1}}$)	(${\mathit{\lambda }}_{\mathrm{2}}$)
	1	414.209 nm	415.535 nm
A	2	435.689 nm	437.015 nm
	3	439.932 nm	441.258 nm
	1	417.126 nm	418.452 nm
B	2	433.037 nm	434.363 nm
	3	442.849 nm	444.175 nm

Figure 2

Example of spectral radiance sensitivity to NO${}_{\mathrm{2}}$ VCD under reference HIS observation conditions at Pohang area shown with three different wavelength pairs per each wavelength pair type (Type_A in red and Type_B in blue) used in the MWP method.

[Figure omitted. See PDF]

Multiple wavelength pairs are selected within the fitting window to reduce the uncertainty. However, when the NO${}_{\mathrm{2}}$ VCDs are retrieved based on the wavelength pair fitting method, which has been applied to the ground-based HIS observations (Park et al., 2019), wavelength pairs of Type_A or Type_B can attribute systematic biases in the retrieved NO${}_{\mathrm{2}}$ VCDs as the assumed surface and aerosol reflectance spectral dependency in the RTM simulation differs from the real world. For example, when the wavelength pairs are only comprised of Type_A, the $R$ value for the wavelength pair is expected to be negatively correlated with the NO${}_{\mathrm{2}}$ VCD. In conditions where the reflectance is higher at the longer wavelength, which is common for land or vegetated surface, the $R$ value is likely to be underestimated, resulting in a positively biased NO${}_{\mathrm{2}}$ VCD. On the other hand, conditions with lower reflectance at a longer wavelength, which is a general feature over water, lead to an overestimation of the $R$ value and an underestimation of the NO${}_{\mathrm{2}}$ VCD. To mitigate these biases, it is necessary to consider accurate and highly resolved spectral surface reflectance and the aerosol optical property data, temporally and spatially collocated to the research flights. However, acquiring such data is often challenging, and even with massive computational consumption (i.e., running high-resolution CTMs), it is likely to be uncertain.

Given the limitations in obtaining suitable surface and aerosol data, a new analytical solution has been derived aiming to minimize the influence of the spectral dependency of surface reflectivity and aerosol optical properties. For a single analytical solution (i.e., NO${}_{\mathrm{2}}$ VCD), a set of wavelength pairs comprised of one Type_A and one Type_B wavelength pair was used to further reduce the influence of the random uncertainty. Therefore, three independent analytic solutions from three sets of wavelength pairs yield three NO${}_{\mathrm{2}}$ VCD values from a single observed spectrum. Specifically, a set of wavelength pairs was formed by grouping the neighboring wavelength pairs with the same index (shown in Table 2). 1$$\left\{\begin{array}{l}{\mathrm{\text{VCD}}}_{\mathrm{A}}=R_{\mathrm{rtm},\mathrm{A}}\cdot a_{\mathrm{A}}+b_{\mathrm{A}},a_{\mathrm{A}}<\mathrm{0}\left(\text{VCD}\uparrow \to R\downarrow \right)\\ {\mathrm{\text{VCD}}}_{\mathrm{B}}=R_{\mathrm{rtm},\mathrm{B}}\cdot a_{\mathrm{B}}+b_{\mathrm{B}},a_{\mathrm{B}}>\mathrm{0}\left(\text{VCD}\uparrow \to R\uparrow \right)\end{array}\right.$$ Equation (1) shows the relationship between NO${}_{\mathrm{2}}$ VCD (abbreviated as “VCD” in the equations for simplicity) and the $R$ values of both Type_A and Type_B wavelength pairs, assuming a first-order linear relationship between the two. Here, $a_{\mathrm{A}}$ and $a_{\mathrm{B}}$ represent the slope, while $b_{\mathrm{A}}$ and $b_{\mathrm{B}}$ represent the intercept. This assumption was confirmed from the forward RTM simulations (see details in Sect. 3.2) by the squared correlation coefficient ($r^{\mathrm{2}}$) of the regressed line (i.e., Eq. 1) exceeding 0.99 (shown in Fig. S1 in the Supplement) for all six wavelength pairs.

As we assume the ratio between real-world reflectivity (i.e., from the surface and aerosols) and the assumed reflectivity for the RTM at ${\mathit{\lambda }}_{\mathrm{2}}$ is a $k$ fold of the value at ${\mathit{\lambda }}_{\mathrm{1}}$, the relation between the $R$ value of the wavelength pair simulated from the RTM and the expected $R$ value calculated from real-world HIS observation can be expressed as Eq. (2). The same $k$ value can be assumed for the wavelength pairs grouped as the same wavelength pair set, since the difference between the wavelengths of each wavelength pair ($\mathrm{\Delta }\mathit{\lambda }={\mathit{\lambda }}_{\mathrm{2}}-{\mathit{\lambda }}_{\mathrm{1}}$) is small ($<\mathrm{1.5}$ nm) and almost consistent at all six wavelength pairs. Moreover, the wavelengths constitute a single wavelength pair set range within a 5 nm bandwidth, so we can assume the spectral dependency to be invariant. 2$$\left\{\begin{array}{l}{R}_{\mathrm{obs},\mathrm{A}}=\frac{\mathrm{1}}{k}\cdot R_{\mathrm{rtm},\mathrm{A}}\\ R_{\mathrm{obs},\mathrm{B}}={\displaystyle \frac{\mathrm{1}}{k}}\cdot R_{\mathrm{rtm},\mathrm{B}}\end{array}\right.$$ Here, subscripts “rtm” and “obs” refer to the simulated value from the RTM and the expected value from the real hyperspectral observations, respectively. From Eqs. (1) and (2), the relation between the biased VCD estimates (VCD${}_{\mathrm{obs}}$) and the observed radiance ratio of the wavelength pair ($R_{\mathrm{obs}}$) can also be expressed with the simulated radiance ratio ($R_{\mathrm{rtm}}$) together with the spectral dependency factor $k$ (Eq. 3). 3$$\left\{\begin{array}{l}{\text{VCD}}_{\mathrm{obs},\mathrm{A}}=R_{\mathrm{obs},\mathrm{A}}\cdot a_{\mathrm{A}}+b_{\mathrm{A}}=R_{\mathrm{rtm},\mathrm{A}}\cdot \frac{a_{\mathrm{A}}}{k}+b_{\mathrm{A}}\\ {\text{VCD}}_{\mathrm{obs},\mathrm{B}}=R_{\mathrm{obs},\mathrm{B}}\cdot a_{\mathrm{B}}+b_{\mathrm{B}}=R_{\mathrm{rtm},\mathrm{B}}\cdot {\displaystyle \frac{a_{\mathrm{B}}}{k}}+b_{\mathrm{B}}\end{array}\right.$$ Since the simulated NO${}_{\mathrm{2}}$ VCDs should be the same regardless of wavelength pair type, the “true” NO${}_{\mathrm{2}}$ VCD (VCD${}_{\mathrm{True}}$) can be expressed as follows (Eq. 4). 4$$\begin{array}{rl}{\text{VCD}}_{\mathrm{True}}& ={\text{VCD}}_{\mathrm{rtm},\mathrm{A}}={\text{VCD}}_{\mathrm{rtm},\mathrm{B}}=R_{\mathrm{rtm},\mathrm{A}}\cdot a_{\mathrm{A}}+b_{\mathrm{A}}\\ & =R_{\mathrm{rtm},\mathrm{B}}\cdot a_{\mathrm{B}}+b_{\mathrm{B}}\end{array}$$ From Eq. (4), $R_{\mathrm{rtm},\mathrm{A}}$ can be shown as Eq. (5). 5$$R_{\mathrm{rtm},\mathrm{A}}={\displaystyle \frac{R_{\mathrm{rtm},\mathrm{B}}\cdot a_{\mathrm{B}}+b_{\mathrm{B}}-b_{\mathrm{A}}}{a_{\mathrm{A}}}}$$ Then the first equation of the simultaneous equations of Eq. (3) can be reformulated by replacing $R_{\mathrm{rtm},\mathrm{A}}$ with Eq. (5), yielding Eq. (6). 6$${\text{VCD}}_{\mathrm{obs},\mathrm{A}}-b_{\mathrm{A}}={\displaystyle \frac{R_{\mathrm{rtm},\mathrm{B}}\cdot a_{\mathrm{B}}+b_{\mathrm{B}}-b_{\mathrm{A}}}{k}}$$ Equations (3) and (6) can be reorganized to Eq. (7) and then be solved to a single equation (Eq. 8). 7$$\left\{\begin{array}{l}\frac{\mathrm{1}}{k}=\frac{{\text{VCD}}_{\mathrm{obs},\mathrm{A}}-b_{\mathrm{A}}}{R_{\mathrm{rtm},\mathrm{B}}\cdot a_{\mathrm{B}}+b_{\mathrm{B}}-b_{\mathrm{A}}}\\ {\displaystyle \frac{\mathrm{1}}{k}}={\displaystyle \frac{{\text{VCD}}_{\mathrm{obs},\mathrm{B}}-b_{\mathrm{B}}}{R_{\mathrm{rtm},\mathrm{B}}\cdot a_{\mathrm{B}}}}\end{array}\right.$$ 8$${\displaystyle \frac{{\text{VCD}}_{\mathrm{obs},\mathrm{A}}-b_{\mathrm{A}}}{R_{\mathrm{rtm},\mathrm{B}}\cdot a_{\mathrm{B}}+b_{\mathrm{B}}-b_{\mathrm{A}}}}={\displaystyle \frac{{\text{VCD}}_{\mathrm{obs},\mathrm{B}}-b_{\mathrm{B}}}{R_{\mathrm{rtm},\mathrm{B}}\cdot a_{\mathrm{B}}}}$$ Equation (8) can be reorganized as an expression of $R_{\mathrm{rtm},\mathrm{B}}$, yielding Eq. (9). 9$$\begin{array}{rl}{R}_{\mathrm{rtm},\mathrm{B}}& ={\displaystyle \frac{\mathrm{1}}{a_{\mathrm{B}}}}\cdot {\displaystyle \frac{b_{\mathrm{A}}-b_{\mathrm{B}}}{{\text{VCD}}_{\mathrm{obs},\mathrm{A}}-b_{\mathrm{A}}}}\\ & \cdot {\displaystyle \frac{\mathrm{1}}{\left(\frac{\mathrm{1}}{{\text{VCD}}_{\mathrm{obs},\mathrm{A}}-b_{\mathrm{A}}}-\frac{\mathrm{1}}{{\text{VCD}}_{\mathrm{obs},\mathrm{B}}-b_{\mathrm{B}}}\right)}}\\ & ={\displaystyle \frac{\mathrm{1}}{a_{\mathrm{B}}}}\cdot {\displaystyle \frac{\left(b_{\mathrm{A}}-b_{\mathrm{B}}\right)\cdot \left({\text{VCD}}_{\mathrm{obs},\mathrm{B}}-b_{\mathrm{B}}\right)}{{\text{VCD}}_{\mathrm{obs},\mathrm{B}}-{\text{VCD}}_{\mathrm{obs},\mathrm{A}}+b_{\mathrm{A}}-b_{\mathrm{B}}}}\end{array}$$ Then the final equation, Eq. (10), for the “true” unbiased NO${}_{\mathrm{2}}$ VCD is derived only as a function of the observed $R$ values of the Type_A and Type_B wavelength pairs ($R_{\mathrm{obs},\mathrm{A}}$, $R_{\mathrm{obs},\mathrm{B}}$) and the coefficients for the first-order polynomials regressed upon the relationship between simulated $R$ values and the NO${}_{\mathrm{2}}$ VCDs (i.e., $a_{\mathrm{A}}$, $b_{\mathrm{A}}$, $a_{\mathrm{B}}$, and $b_{\mathrm{B}}$ in Eq. 1). 10$$\begin{array}{rl}{\text{VCD}}_{\mathrm{True}}& =R_{\mathrm{rtm},\mathrm{B}}\cdot a_{\mathrm{B}}+b_{\mathrm{B}}\\ & ={\displaystyle \frac{\left(b_{\mathrm{A}}-b_{\mathrm{B}}\right)\cdot \left({\text{VCD}}_{\mathrm{obs},\mathrm{B}}-b_{\mathrm{B}}\right)}{{\text{VCD}}_{\mathrm{obs},\mathrm{B}}-{\text{VCD}}_{\mathrm{obs},\mathrm{A}}+b_{\mathrm{A}}-b_{\mathrm{B}}}}+b_{\mathrm{B}}\\ & ={\displaystyle \frac{\left(b_{\mathrm{A}}-b_{\mathrm{B}}\right)\cdot \left(R_{\mathrm{obs},\mathrm{B}}\cdot a_{\mathrm{B}}\right)}{R_{\mathrm{obs},\mathrm{B}}\cdot a_{\mathrm{B}}-R_{\mathrm{obs},\mathrm{A}}\cdot a_{\mathrm{A}}}}+b_{\mathrm{B}}\end{array}$$ As shown in Eq. (10), in which the NO${}_{\mathrm{2}}$ VCD is expressed without the broadband spectral dependency term ($k$), the NO${}_{\mathrm{2}}$ VCDs can be retrieved with the MWP method without considering the varying spectral variability in the surface or the aerosol properties. For the MWP method application, a set of wavelength pairs, one Type_A and one Type_B, is required. Therefore, from three sets of wavelength pairs, three different NO${}_{\mathrm{2}}$ VCD values can be retrieved from a single observed spectrum. It is noteworthy to mention that the spectral binning was applied to the HIS-observed spectra and subsequently the RTM-simulated radiance spectra before implementing the MWP method, aiming to increase the signal-to-noise ratio (SNR). The spectral binning was constrained within a range that does not sacrifice the sensitivity of the wavelength pair method (Fig. 2) and takes into account the HIS oversampling rate, resulting in spectral binning of $\pm \mathrm{2}$ original spectral bins.

3.2 Forward radiative transfer modeling

In order to apply the MWP method to the HIS-observed spectra, a look-up table (LUT) of the $R$ values corresponding to different NO${}_{\mathrm{2}}$ VCDs under various environments and observation conditions should be constructed from the forward RTM simulations for each wavelength pair. UVSPEC, specializing in radiative transfer calculations in UV–VIS bands, within the versatile libRadtran (version 2.0.4) radiative transfer software package (Emde et al., 2016; Mayer and Kylling, 2005) was utilized for the forward radiative transfer simulations. Spectral radiances were simulated from UVSPEC based on input meteorological parameters and the profiles of various atmospheric constituents such as trace gases and aerosols.

Table 3 shows the input conditions for UVSPEC for the LUT construction. Some of the input entries were adjusted to each research flight to account for possible conditions that vary in location and time. The ECMWF Reanalysis v5 (ERA5) hourly reanalysis data (Hersbach et al., 2020, 2023a) for the corresponding flight area and time of each research flight were used as an input atmospheric condition (i.e., vertical profiles of pressure, temperature, and mixing ratio of major atmospheric gaseous constituents except for NO${}_{\mathrm{2}}$) for the RTM simulations. The vertical profile of NO${}_{\mathrm{2}}$ can significantly alter the AMF or the effective absorption of NO${}_{\mathrm{2}}$ per unit VCD, so it is necessary to use as realistic a profile as possible. The most intuitive way is to use vertical profiles simulated from CTM with sufficiently high spatial resolution. However, unlike conventional AMF calculations using LUT of scattering weights and taking the inner product with CTM-simulated vertical profiles, it is not feasible to take all the vertical profiles within the target domain into account and construct the LUT of the simulated radiances, which the MWP method requires.

Table 3

Variables and corresponding entries for the radiative transfer model simulations to construct the LUT.

Variables (acronym)	Entry range (number of entries)	Units	Remarks
NO${}_{\mathrm{2}}$ VCD	0.15–3.60 (70)	Dobson unit (DU)	1 DU $=\mathrm{2.687}\times {\mathrm{10}}^{\mathrm{16}}$ molec. cm${}^{-\mathrm{2}}$
Planetary boundary layer height (PBLH)	0.2–1.8 (5)	km	–
Surface reflectivity (ALB)	0.01, 0.03, 0.05, 0.08, 0.11, 0.15,	–	Assumed as Lambertian surface
	0.20, 0.25, 0.30, 0.35 (10)		without spectral dependency
Solar zenith angle (SZA)	Variable per flight (8–10)	Degree (${}^{\circ }$)	3${}^{\circ }$ interval, ranging from 42 to 75${}^{\circ }$
			(SZA $>$ 70${}^{\circ }$ are filtered out)
Flight altitude (ALT)	Variable per flight (4–9)	m (above the mean	100 m interval near
		sea level)	the cruising altitude
Viewing zenith angle (VZA)	0–24 (9)	Degree (${}^{\circ }$)	-
Relative azimuth angle (RAA)	0–180 (7)	Degree (${}^{\circ }$)	–
Wavelength	400.15–459.82 (226)	nm	Encompassing the fitting window

As an alternative, NO${}_{\mathrm{2}}$ vertical profiles from the high-resolution (3 km $\times$ 3 km horizontal resolution; 23 vertical grids ranging from the surface to 70 hPa) Community Multiscale Air Quality (CMAQ) v5.2 model (Appel et al., 2017; US EPA, 2017) simulations at each target domain were collected to calculate representative NO${}_{\mathrm{2}}$ vertical profiles under certain planetary boundary layer height (PBLH) and NO${}_{\mathrm{2}}$ VCD conditions. The GMAP/SIJAQ (Geostationary Environmental Monitoring Spectrometer Satellite Integrated Joint Monitoring of Air Quality) v2.0 emission inventory based on the Clean Air Policy Support System (CAPSS) 2018 and Comprehensive Regional Emissions inventory for Atmospheric Transport Experiment (CREATE) v5.0 (Woo et al., 2020) under the framework of the SMOKE-Asia (Woo et al., 2012) emission model was used for high-resolution (3 km $\times$ 3 km) anthropogenic emissions for the CMAQ, while the Model of Emissions of Gases and Aerosols from Nature (MEGAN) v2.10 was used for natural emissions. Meteorological parameters were simulated from the high-resolution (3 km $\times$ 3 km horizontal resolution; 31 vertical grids ranging from the surface to 50 hPa) Weather Research and Forecasting (WRF) v3.9.1 model (Skamarock et al., 2008), with the initial conditions from the National Centers for Environmental Prediction Final Analysis (NCEP/FNL) reanalysis data. More specifically, NO${}_{\mathrm{2}}$ vertical profiles within each target domain were grouped by the CMAQ-driven PBLH with a $\pm \mathrm{0.2}$ km range from the PBLH entries shown in Table 3. The variability (i.e., quantiles) of NO${}_{\mathrm{2}}$ mixing ratios at each altitude level was calculated from vertical profiles grouped under the same PBLH conditions, and the NO${}_{\mathrm{2}}$ vertical profile corresponding to certain a NO${}_{\mathrm{2}}$ VCD was deduced from Levenberg–Marquardt fitting with the quantile value that best fits the targeted column-integrated NO${}_{\mathrm{2}}$ VCD (Fig. S2). It should be noted that the stratospheric NO${}_{\mathrm{2}}$ VCD (i.e., NO${}_{\mathrm{2}}$ column densities above 15 km in altitude) was assumed constant (0.185 DU) based on the US standard atmosphere (1976) since the CMAQ only spans altitudes up to approximately 70 hPa.

Appropriate consideration of the real-world aerosol vertical profile is as important as accurately depicting the NO${}_{\mathrm{2}}$ vertical profile because the aerosol distribution compels the effective light path (or the scattering weight). However, due to the lack of a vertically resolved observation dataset, column-integrated aerosol optical properties such as aerosol optical depth (AOD) and single scattering albedo (SSA) from the nearest AERONET (Holben et al., 1998; Kim et al., 2007; Table S1 in the Supplement) sites were used for the radiative transfer calculation.

For the forward radiance simulation with UVSPEC, the TSIS-1 hybrid solar reference spectrum (Bak et al., 2021; Coddington et al., 2021) was used for the solar irradiance spectrum, while the surface was assumed to exhibit Lambertian reflectance without spectral dependency. Input spectra such as absorption cross-sections of trace gases and the solar irradiance spectrum were convolved by using the Gaussian slit function with FWHM of 1.4 nm to ensure conformity with the expected observed spectra from the HIS.

3.3 Spectral and radiometric calibrations

Accurate spectral and radiometric characterization of an instrument is essential to assure consistency and minimize the uncertainty of the retrieved NO${}_{\mathrm{2}}$ VCDs. Therefore, the spectral and radiometric properties of the HIS have been assessed, as well as post-calibrations upon the obtained spectra from the airborne HIS observations. First, the SNR of the HIS was examined based on dark samples (i.e., HIS-observed spectra without light exposure) collected with diverse exposure times. Dark signals (CCD signals from the dark samples) are composed of dark current – a thermally induced false signal – and the offset value – an arbitrary value added by the manufacturer to avoid unintended positive bias from instrument noise. Dark signals had almost no dependency on the exposure time, implying minimal dark currents of the HIS, and the dark offset ranged from 380 to 410 counts from spatial column to column. However, the noise level of the pixels of interest (spectral rows within the fitting window used in the MWP method) was substantial, necessitating additional binning to increase the SNR by controlling the noise. Since the optimal spectral binning had already been determined and addressed in the MWP algorithm (discussed in Sect. 3.1), 50 consecutive CCD spatial columns were binned to a single co-added column, starting from the centerline column, yielding 27 co-added spatial columns and discarding 21 raw lateral columns on each side. The noise level estimated after the spatial binning was approximately 2–2.5 counts under an integration time of 1 s (Fig. S3), which would result in an SNR ranging on the scale of a few thousand (depending on the signal intensity) during the airborne observations. Spatial binning was implemented in the HIS-observed spectra and has been taken into account for all the subsequent calibrations below.

The spectral calibrations of the instrument and the observed spectra from the HIS include (1) slit function characterization, (2) initial spectral registration for spatially binned CCD pixels, and (3) the retrieval of wavelength shifts from the initial spectral registration for every observed frame and spatial column. The slit function characterization and initial spectral registration were carried out by illuminating the HIS with a reference Hg lamp, which emits pseudo-monochromatic emission lines from the low-pressure Hg gas in the UV–VIS range. The shape of the slit function was assumed to be a simple Gaussian distribution, and the slit function FWHM was retrieved for each Hg-lamp emission line within the HIS spectral range. The slit function FWHM generally conformed to the instrument specifications (Table 1) but varied within the range of 1.1–1.4 nm. The variability in the slit function FWHM primarily arises from the deviation of the actual slit function from the Gaussian distribution due to CCD non-uniformity and the instability in the optics (slit, grating, etc.) caused by thermal variation and vibrations. The initial spectral registration of the spatially binned CCD pixels was also accomplished by first finding the spectral row numbers corresponding to the centerline of the illuminated CCD pixels for each Hg-lamp emission line with known wavelengths. Then, the wavelengths for the rest of the spectral rows were registered by interpolating the spectral rows with known wavelengths.

The initial pre-flight wavelength registration was calibrated since the HIS is exposed to varying temperatures during the flight, which can result in a spectral shift from the initial spectral registration. The average spectral shift within the fitting window (410–460 nm) of the MWP method was calculated for every HIS-observed frame and spatial columns by fitting the spatially binned HIS spectra with the convolved solar irradiance spectrum using Levenberg–Marquardt nonlinear fitting. Polynomials for scaling and offset terms, together with spectral shift and the FWHM value of the slit function, were considered for the fitting. The detailed descriptions of spectral calibration, including slit function characterization, initial spectral registration, and spectral shift calculation, can be found in previous papers (Kang et al., 2020; Liu et al., 2015; Nowlan et al., 2016).

In order to convert raw CCD counts into a radiance, it is necessary to obtain precise pixel-by-pixel sensitivity from specialized facilities equipped with a well-calibrated reference lamp and integrating sphere. However, the intension of using the low-cost HIS is to make the calibration and maintenance of the instrument as simple as possible for feasible observations, and equipping them in such costly facilities would contradict this. Moreover, the radiometric sensitivity of the CCD is altered through time and with temperature change, and even the most accurate absolute radiometric calibration does not guarantee the integrity of the HIS-observed radiance during the flight. Therefore, an alternative approach based on HIS spectra obtained from the clean pixel was used to minimize the effect of non-uniformity in CCD sensitivities.

The basic concept of clean-pixel calibration is determining the actual in-flight pixel-by-pixel CCD sensitivities by comparing the HIS spectra (i.e., raw digital counts) at the clean pixel with the simulated radiance from the RTM. The spectral scale factor (SF), a spectral conversion factor corresponding to the CCD sensitivity, is calculated for every pre-binned spatial column by dividing the RTM-simulated radiance by the HIS-observed raw counts. It is important to accurately account for the parameters affecting the effective light absorption of NO${}_{\mathrm{2}}$, such as surface reflectivity and NO${}_{\mathrm{2}}$ VCD, in the RTM simulation over the clean pixel. Therefore, clean pixels were selected as a place with easily identifiable homogeneous surface conditions (i.e., ocean) and with a background level of NO${}_{\mathrm{2}}$ near the surface. It is also important to select a clean pixel collected after a sufficient time for instrument stabilization at the observation altitude to avoid sensitivity change due to temperature change.

Applying spectral SF as an alternative radiometric calibration can result in a substantial error in the absolute magnitude of converted radiances. However, the MWP method uses the relative radiance ratio of the wavelength pair, and the pixel-to-pixel non-uniformity is the only primary concern, not the absolute magnitude of the conversion factor. Moreover, by applying SF calculated as a ratio between simulated and observed spectra, the discrepancies between the RTM simulation settings and the real world, including the uncertainties in instrument characterization, can be corrected empirically. To take this advantage, it is necessary to retrieve spectral SF after accurate spectral calibration and to use respective SF for every research flight, different flying altitude, and exposure time. Meanwhile, on the other side of the virtue of applying the spectral SF from a clean pixel, it has one critical drawback in estimating the NO${}_{\mathrm{2}}$ VCD value corresponding to the clean pixel. Whenever a certain NO${}_{\mathrm{2}}$ VCD value is determined for a clean pixel, it directly affects all the other NO${}_{\mathrm{2}}$ VCDs retrieved from the HIS using the spectral SF obtained with a particular clean pixel. As a result, a radiometric calibration using spectral SF yields relative NO${}_{\mathrm{2}}$ VCD, which may be biased proportionally to the bias of NO${}_{\mathrm{2}}$ VCD estimates at the clean pixel. Nevertheless, the application of spectral SF offers more significant merits than downsides, especially for calibrating low-cost hyperspectral sensors with unstable optical-radiometric characteristics that are operated without sophisticated laboratory calibrations or maintenance since it not only calibrates pixel non-uniformity but also addresses other uncertainties arising from the inaccurate depiction of instrumental characteristics during the flight.

3.4 Uncertainty estimation

Uncertainties of the NO${}_{\mathrm{2}}$ VCDs retrieved from the MWP method can be estimated from the errors in input parameters and their sensitivity to the resulting NO${}_{\mathrm{2}}$ VCD. One of the advantages of the MWP method is that the solution (NO${}_{\mathrm{2}}$ VCD) is retrieved analytically, making the uncertainty estimation feasible by adopting Gaussian error propagation. The error estimation of the NO${}_{\mathrm{2}}$ VCDs starts from the final equation (Eq. 10) of the MWP method discussed in Sect. 3.1.

The regression parameters for the wavelength pair set (i.e., $a_{\mathrm{A}}$, $b_{\mathrm{A}}$, $a_{\mathrm{B}}$, and $b_{\mathrm{B}}$ in Eq. 10) are calculated from the NO${}_{\mathrm{2}}$ VCD to $R$ value relationship, whereas the $R$ values for the regression are retrieved by interpolating the LUT with the corresponding observation conditions (i.e., SZA, VZA, RAA, PBLH). Therefore, these regression parameters can be assumed as assured values for a certain observation condition. However, when the real-world observation conditions differ from the input atmospheric conditions assumed for the retrieval, the $R$ values will change, posing an error in the resulting NO${}_{\mathrm{2}}$ VCD. Accordingly, the sensitivity of these $R$ values as a response to the uncertainties of input parameters was calculated in order to estimate the error propagated to the final NO${}_{\mathrm{2}}$ VCD retrieval in the MWP method.

To calculate error propagation, $R_{\mathrm{obs}}$ values in Eq. (10) should be transposed with $R_{\mathrm{rtm}}$, assuming the spectral dependency term $k$ in Eq. (2) holds constant within the wavelength pair set. Then the pre-constructed LUT of $R$ values can be utilized to examine the sensitivity of $R$ values to errors in input parameters and can be used in estimating the uncertainty range of the NO${}_{\mathrm{2}}$ VCD. Furthermore, we can transpose the ratio between $R_{\mathrm{rtm},\mathrm{A}}$ and $R_{\mathrm{rtm},\mathrm{B}}$ as $Q$ and reformulate Eq. (10) as the univariate equation of $Q$ (Eq. 11). This transposition is necessary to avoid the necessity of examining the correlation between $R_{\mathrm{rtm},\mathrm{A}}$ and $R_{\mathrm{rtm},\mathrm{B}}$ for the error propagation.

11$$\begin{array}{rl}{\text{VCD}}_{\mathrm{True}}& ={\displaystyle \frac{a_{\mathrm{B}}\cdot \left(b_{\mathrm{A}}-b_{\mathrm{B}}\right)}{a_{\mathrm{B}}-\left(\frac{R_{\mathrm{obs},\mathrm{A}}}{R_{\mathrm{obs},\mathrm{B}}}\right)\cdot a_{\mathrm{A}}}}+b_{\mathrm{B}}\\ & ={\displaystyle \frac{a_{\mathrm{B}}\cdot \left(b_{\mathrm{A}}-b_{\mathrm{B}}\right)}{a_{\mathrm{B}}-\left(\frac{R_{\mathrm{rtm},\mathrm{A}}}{R_{\mathrm{rtm},\mathrm{B}}}\right)\cdot a_{\mathrm{A}}}}+b_{\mathrm{B}}\\ & ={\displaystyle \frac{a_{\mathrm{B}}\cdot \left(b_{\mathrm{A}}-b_{\mathrm{B}}\right)}{a_{\mathrm{B}}-Q\cdot a_{\mathrm{A}}}}+b_{\mathrm{B}}\end{array}$$ Under these circumstances, the error in VCD${}_{\mathrm{True}}$ can be estimated from Eq. (11) as follows (Eq. 12). 12$${\mathit{\sigma }}_{{\mathrm{VCD}}_{\mathrm{True}}}^{\mathrm{2}}=\left({\displaystyle \frac{a_{\mathrm{A}}\cdot a_{\mathrm{B}}\cdot \left(b_{\mathrm{A}}-b_{\mathrm{B}}\right)}{{\left(a_{\mathrm{B}}-\cdot a_{\mathrm{A}}\right)}^{\mathrm{2}}}}\right)\cdot {\mathit{\sigma }}_Q^{\mathrm{2}}$$ The uncertainty of the $Q$ value (${\mathit{\sigma }}_Q$) is determined by the sensitivity of $Q$ values to the uncertainty of input parameters. However, ${\mathit{\sigma }}_Q$ has a dependency on $Q$ values, so the uncertainty of the normalized $Q$ value (${\mathit{\sigma }}_{Q_{\mathrm{rel}}}$), which better portrays the relative uncertainty, has been calculated from Eq. (13) under the assumption that all the parameters are uncorrelated. It is noteworthy to mention that the ${\mathit{\sigma }}_{Q_{\mathrm{rel}}}$ can be easily converted to ${\mathit{\sigma }}_Q$ by multiplying the $Q$ value calculated from the HIS spectrum in the actual NO${}_{\mathrm{2}}$ VCD retrieval. 13$${\mathit{\sigma }}_{Q_{\mathrm{rel}}}^{\mathrm{2}}={\sum }_{\mathrm{var}}\left({\left|{\displaystyle \frac{\partial Q_{\mathrm{norm}}}{\partial \text{var}}}\right|}^{\mathrm{2}}\cdot {\mathit{\sigma }}_{\mathrm{var}}^{\mathrm{2}}\right)$$ Table 4 presents the assumed uncertainty range of input parameters (${\mathit{\sigma }}_{\mathrm{var}}$) used in estimating ${\mathit{\sigma }}_{Q_{\mathrm{rel}}}$, and Fig. 3 (also, see Figs. S4 and S5) is an example of the sensitivity test results showing $\frac{\partial Q_{\mathrm{norm}}}{\partial \text{var}}{\mathit{\sigma }}_{\mathrm{var}}$ values under the reference conditions (presented in Table S2). Most of the uncertainties were posed by the uncertainties in the calibrations, with the greatest contribution from the uncertainties of spectral SF calibration followed by the uncertainties in spectral shift calibration. The uncertainty of the spectral SF calibration largely depended on the quality of the clean pixel used for the spectral SF derivation. For instance, when the spectral SF is calculated with clean-pixel data obtained over a homogeneous surface and with the nearly steady-signal state, the uncertainty propagated to the $Q$ value is minimal. The general relative uncertainty of a spectral SF in converting the raw HIS signal to the non-uniformity-corrected spectra was kept below ${\mathrm{10}}^{-\mathrm{2}}$ and ${\mathrm{10}}^{-\mathrm{1}}$ for the clean pixel over the ocean and the land surface, respectively. The accuracy of the spectral calibration was generally high, with the spectral shift uncertainty constrained within 0.05 nm on most of the occasions. However, $Q$ values were even sensitive to the small spectral shifts, resulting in the second largest contributor to the total uncertainty.

Table 4

Uncertainty range of variables (${\mathit{\sigma }}_{\mathrm{var}}$) assumed for the error estimation.

Variables	Uncertainty	Unit	Remarks
	(${\mathit{\sigma }}_{\mathrm{var}}$)
PBLH	0.4	km	–
Albedo (reflectivity)	0.03	–	–
Solar zenith angle	0.05	Degree (${}^{\circ }$)	–
Flight altitude	10	m	–
Viewing zenith angle	3	Degree (${}^{\circ }$)	–
Relative azimuth angle	3	Degree (${}^{\circ }$)	Uncertainty mainly attributed to viewing azimuth angle
Wavelength shift	0.05	nm	Wavelength shift calibration fitting uncertainty
Instrumental noise	Variable	–	Depending on the observed signal intensity
Radiometric calibration	Variable	–	Depending on which spectral SF applied for the retrieval

Figure 3

Sensitivity of $Q$ values calculated from simulated radiances at wavelength pair 1 (i.e., Type_A: 414.209, 415.535 nm; Type_B: 417.126, 418.452 nm) depending on (a) NO${}_{\mathrm{2}}$ VCD, (b) PBLH, (c) albedo (ALB; reflectivity), (d) solar zenith angle (SZA), (e) observation altitude (ALT), (f) viewing zenith angle (VZA), and (g) relative azimuth angle (RAA) considering atmospheric condition at Pohang and its corresponding reference conditions (shown in Table S2).

[Figure omitted. See PDF]

The relative significance of the instrument noise level is characterized by SNR and is dependent on the signal intensity measured by the HIS. Assuming an SNR of $\mathrm{2}\times {\mathrm{10}}^{\mathrm{3}}$, which is a typical value in between the scenes over the bright land surface and the dark ocean surface, caused the instrument noise to become a subsequent source of error.

Among the variables composing the LUT for the MWP method, surface reflectance and PBLH were the most influential parameters in terms of uncertainties. PBLH from ERA5 reanalysis data (Hersbach et al., 2023b) was adopted for an application of the MWP method because there was only a subset of research flights having explicit observations of PBLH, such as from lidar, located in the vicinity of the target domain. The PBLHs from lidar observations and the ERA5 reanalysis were comparable to those from the research flights with the lidar observation nearby (not shown), but the uncertainty of region-specific PBLH in the global reanalysis dataset was regarded conservatively considering that the target domains are mostly located near the coastline where the PBLH can exhibit significant spatial variability.

The Suomi NPP VIIRS BRDF data (Suomi National Polar-orbiting Partnership Visible Infrared Imaging Radiometer Suite Bidirectional Reflectance Distribution Function; VNP43C1; Schaaf et al., 2019) was used to calculate effective surface reflectivity accounting for the viewing geometry during the airborne HIS observation. Moreover, HIS observations with extreme viewing zenith angle were discarded in the retrieval, leaving only a small space for an error in surface reflectance. However, the spatial inhomogeneity of surface reflectance can be significant; hence the uncertainty of the surface reflectance should be considered. The sensitivity of $Q$ values to the surface reflectance was greater at the darker surface, probably attributed to the perturbation amplified in terms of relative change under low-reflectivity conditions due to a smaller denominator value. It is pertinent to mention that the uncertainties in $Q$ values, induced by the uncertainties in PBLH and surface reflectance, become greater under high NO${}_{\mathrm{2}}$ VCD conditions. NO${}_{\mathrm{2}}$ becomes more concentrated at the low altitudes as the column density increases, coupled with the fact that both PBLH and surface reflectivity control the scattering weights of the low altitudes, resulting in a greater sensitivity of the $Q$ values to the PBLH and surface reflectivity perturbations under high NO${}_{\mathrm{2}}$ VCD conditions.

Solar zenith angle and the observing altitudes are almost certain, considering the accuracy of chronographs and the GPS sensor used for the research flights. Nonetheless, the time synchronization between the instruments might have a small discrepancy, and the narrow window of uncertainty has been accounted for in the sensitivity test. Furthermore, the viewing zenith angle and the relative azimuth angle can show moderate uncertainty from the misalignment between the HIS and the IMU (inertial measurement unit) sensor module, whereas both variables showed a minimal impact on the $Q$ value uncertainty.

The LUT of the sensitivities to uncertainties of each input parameter ($\frac{\partial Q_{\mathrm{norm}}}{\partial \text{var}}{\mathit{\sigma }}_{\mathrm{var}}$) has been constructed for every set of wavelength pairs, incorporating the conditions that the $R$ value LUT spans, for an uncertainty estimation of NO${}_{\mathrm{2}}$ VCDs retrieved from the MWP method. However, it is imperative to acknowledge that there can be additional sources of uncertainty which could further increase the proposed uncertainty range in this study. We have not accounted for uncertainties introduced by assuming aerosol properties or the NO${}_{\mathrm{2}}$ vertical profiles, as well as consequential structural uncertainty caused by varying effective temperature of NO${}_{\mathrm{2}}$ (i.e., temperature dependency of NO${}_{\mathrm{2}}$ absorption; Spinei et al., 2014; Vandaele et al., 2003). Quantifying these uncertainties is challenging due to their large and unpredictable variabilities. For instance, estimating the uncertainty caused by the shape of the NO${}_{\mathrm{2}}$ vertical profile is difficult due to its excessive degree of freedom; previous studies assumed it accounts for approximately 10 % of the relative error (Boersma et al., 2018; Tack et al., 2019).

4

Airborne observations of NO${}_{\mathrm{2}}$ VCDs

4.1 Research flights

Three regions in South Korea were selected as the target domain for the research flight: Chungnam, Jecheon, and Pohang. Each domain has one or more industrial point sources with substantial NO${}_{\mathrm{2}}$ emissions (Choo et al., 2023; Kim et al., 2020; Cleansys, 2023), such as power plants, steel yards, petrochemical complexes, and cement kilns. Figure 4 depicts the composite of the quality-assured (i.e., effective cloud fraction $<\mathrm{0.3}$) TROPOMI NO${}_{\mathrm{2}}$ VCD data (version 2.0.4; ESA and KNMI, 2021) spanning 4 years (from 2019 to 2022) for October and November, which aligns with the time frame of the research flights, with the locations of the target domains as well as the precise position of the point sources within each domain. TROPOMI NO${}_{\mathrm{2}}$ VCDs generally displayed similar concentrations compared to the year-round average (Park et al., 2022), with the highest concentration observed at the southeastern side of the Seoul Metropolitan Area (SMA; 0.85 DU). Although the greatest peak is found over the SMA, noticeable increases in NO${}_{\mathrm{2}}$ VCDs were observed in all three target domains. The NO${}_{\mathrm{2}}$ peaks become even more conspicuous and localized near the major industrial sources as the TROPOMI composite is downscaled to a finer resolution (Fig. 4), providing solid evidence of significant NO${}_{\mathrm{2}}$ emissions from these specific point sources.

Figure 4

Location of major industrial point sources (refer to the legend in the upper left panel for their respective acronyms and symbols) in South Korea targeted for airborne HIS observations overlaid by average TROPOMI NO${}_{\mathrm{2}}$ VCD${}_{\mathrm{Total}}$ in October and November from 2019 to 2022. TROPOMI swath datasets were downscaled to grids of 0.25${}^{\circ }$ $\times$ 0.25${}^{\circ }$ for South Korea and 0.01${}^{\circ }$ $\times$ 0.01${}^{\circ }$ for the target domains (i.e., Chungnam, Jecheon, and Pohang) by calculating average NO${}_{\mathrm{2}}$ VCDs from nearby TROPOMI pixels with the weights inversely proportional to the distance. Blue circles show the location of local cities at each target domain.

[Figure omitted. See PDF]

The details of the research flights for airborne HIS observations and the configurations of the HIS mounted on the airplane are shown in Table S3 and Fig. S6. A total of five research flights were conducted under generally cloud-free conditions, with three flights extensively dedicated to the Chungnam domain. The HIS was mounted on gyro-stabilized mounts in the Cessna 208 aircraft, whereas an external canister was used to mount the HIS on the pylon beneath the Beechcraft 1900D aircraft. The gyro-stabilized mount is expected to have dampened the vibrations for the research flights with the Cessna airplane, while the HIS mounted on an external pylon might have been exposed to a moderate level of vibrations during the flights with the Beechcraft 1900D. The position and observation geometry during the airborne HIS observations were obtained from an IMU/GPS module dedicated to the HIS in principle, but the airplane's IMU/GPS information was used for the research flights with the Beechcraft 1900D after calibration with the difference between HIS viewing geometry and the airplane's plane of reference.

The actual flight path and altitude of each research flight are illustrated in Fig. S7. Flight altitude was limited to 11 000 ft ($\sim$3350 m) due to air traffic regulations in South Korea for unpressurized aircraft cabins. The low flight altitude, combined with the narrow FOV of the HIS (approximately 13${}^{\circ }$), resulted in a narrow observation swath ranging between 340 m (at 5000 ft, $\sim$1520 m) and 750 m (at 11 000 ft, $\sim$3350 m).

4.2

NO${}_{\mathrm{2}}$ VCDs in three industrial areas

The MWP method applied to all three sets of wavelength pairs yields three independent NO${}_{\mathrm{2}}$ VCDs with corresponding uncertainty estimates from a single HIS spectrum. To take full advantage of three independent NO${}_{\mathrm{2}}$ VCDs, the weighted mean value was taken as a single representative NO${}_{\mathrm{2}}$ VCD value for the HIS spectrum from a specific spatial column and frame. The weights were determined to be inversely proportional to the squared uncertainty of each NO${}_{\mathrm{2}}$ VCD from three different wavelength pair sets (Eq. 14), which can minimize the random uncertainty theoretically (Xie et al., 2001). 14$${\mathrm{\text{VCD}}}_{\mathrm{final}}={\sum }_{i=\mathrm{1}}^{\mathrm{3}}\left({\mathit{\alpha }}_i{\text{VCD}}_i\right),\mathrm{\text{where}}{\mathit{\alpha }}_i={\displaystyle \frac{\left(\frac{\mathrm{1}}{{\mathit{\sigma }}_i^{\mathrm{2}}}\right)}{{\sum }_j\left(\frac{\mathrm{1}}{{\mathit{\sigma }}_j^{\mathrm{2}}}\right)}}$$ The narrow swath of the airborne HIS observation resulted in a footprint of pre-binned spatial columns having an extremely short across-track length scale (approximately 27 m at 11 000 ft, $\sim$3350 m, and 13 m at 5000 ft, $\sim$1520 m) compared to the along-track length scale (approximately 100 m under an aircraft ground speed of 360 km h${}^{-\mathrm{1}}$ and 1 s of exposure). Moreover, the remaining uncertainties of the NO${}_{\mathrm{2}}$ VCDs from the pre-binned spatial columns were still excessive; thus the additional post-binning was conducted. To be specific, the weighted average was taken over NO${}_{\mathrm{2}}$ VCDs from all 50 pre-binned spatial columns (i.e., all of the pre-binned across-track pixels) and four consecutive frames (integrated exposure time of 4 s) with consideration of uncertainties of each NO${}_{\mathrm{2}}$ VCD (as shown in Eq. 14). Therefore, the final HIS data have a footprint almost symmetric at an altitude of 6000 ft ($\sim$1830 m; 400 m $\times$ 400 m) but can be variable according to the ground speed of an airplane and the flying altitude.

The NO${}_{\mathrm{2}}$ VCDs from the airborne HIS observations and their uncertainties are shown in Figs. 5, 6, and 7, together with the TROPOMI swath data of the same day. The NO${}_{\mathrm{2}}$ VCDs observed from three research flights over the Chungnam area showed different spatial patterns due to different wind fields (Fig. 5). The highest NO${}_{\mathrm{2}}$ concentration was observed on 24 November 2022 when the wind was calm ($<\mathrm{2}$ m s${}^{-\mathrm{1}}$) within the boundary layer (Fig. 5e). Since the atmosphere was stagnant near the surface, the peaks of NO${}_{\mathrm{2}}$ VCDs were observed right above the major emission point sources. The highest concentration was found over the steel yard, with the peak HIS NO${}_{\mathrm{2}}$ VCD reaching up to 2.0 DU, while the same peak from the TROPOMI swath showed 1.0 DU. Likewise, the HIS captured peaks with higher NO${}_{\mathrm{2}}$ VCDs over the petrochemical complex ($>\mathrm{1.2}$ DU), while the NO${}_{\mathrm{2}}$ VCD of the corresponding TROPOMI footprint remained relatively low (0.9 DU). While the small elevated NO${}_{\mathrm{2}}$ VCD signals were discernible over some of the power plants from the airborne HIS observations, TROPOMI failed to capture those minor peaks.

Figure 5

(a) Downscaled average TROPOMI NO${}_{\mathrm{2}}$ VCDs in October and November from 2019 to 2022 and (b) Google Earth image (© Google Earth) of the Chungnam domain with major industrial point sources of NO${}_{\mathrm{2}}$ and populated urban areas denoted with red and blue circles, respectively. Panels (c), (e), and (g) show NO${}_{\mathrm{2}}$ VCD retrieval results from HIS overlaid on TROPOMI swath data, and panels (d), (f), and (h) show the uncertainty of HIS-retrieved NO${}_{\mathrm{2}}$ VCD. Panels (c) and (d), (e) and (f), and (g) and (h) show observed results from the research flight on 17 October 2020, 24 November 2022, and 25 November 2022, respectively.

[Figure omitted. See PDF]

Figure 6

(a) Downscaled average TROPOMI NO${}_{\mathrm{2}}$ VCDs in October and November from 2019 to 2022 and (b) Google Earth image (© Google Earth) of the Jecheon domain with major industrial point source of NO${}_{\mathrm{2}}$ and populated urban area denoted with red and blue circles, respectively. Panels (c) and (d) show observed results from the research flight on 3 November 2020, while panel (c) shows NO${}_{\mathrm{2}}$ VCD retrieval results from HIS overlaid on TROPOMI swath data, and panel (d) shows the uncertainty of HIS-retrieved NO${}_{\mathrm{2}}$ VCD.

[Figure omitted. See PDF]

Figure 7

(a) Downscaled average TROPOMI NO${}_{\mathrm{2}}$ VCDs in October and November from 2019 to 2022 and (b) Google Earth image (© Google Earth) of the Pohang domain with an industrial point source of NO${}_{\mathrm{2}}$ and populated urban area denoted with red and white circles, respectively. Panels (c) and (d) show observed results from the research flight on 5 November 2020, while panel (c) shows NO${}_{\mathrm{2}}$ VCD retrieval results from HIS overlaid on TROPOMI swath data, and panel (d) shows the uncertainty of HIS-retrieved NO${}_{\mathrm{2}}$ VCD.

[Figure omitted. See PDF]

The research flights on 17 October 2020 (Fig. 5c) and 25 November 2022 (Fig. 5g) exhibited relatively lower NO${}_{\mathrm{2}}$ VCDs than the flight on 24 November 2022, attributed to relatively strong surface wind rapidly dissipating the emitted plumes from the point sources. For instance, Fig. 5g shows a clear depiction of an NO${}_{\mathrm{2}}$ plume conically dispersed and transported to the north under clear southerly wind at the boundary layer. The significance of NO${}_{\mathrm{2}}$ plume emitted from a steel yard is also emphasized in Fig. 5g, whereas the stronger wind spreading out the plume into a wider area resulted in lower NO${}_{\mathrm{2}}$ VCDs compared to the previous day (Fig. 5e). Meanwhile, the observation area of the research flight on 17 October 2020 (Fig. 5c) was relatively narrow compared to other flights at Chungnam, only incorporating the petrochemical complex as a major emission point source within the HIS observation area. The wind direction was between northerly and northwesterly, with moderate wind speed at the surface, and the NO${}_{\mathrm{2}}$ plume emitted from the petrochemical complex was clearly distinguished with a dispersing pattern toward the southeast.

It is worth pointing out that in Fig. 5g, TROPOMI shows elevated NO${}_{\mathrm{2}}$ VCDs on the eastern side of the plume detected by the HIS, originating from a steel yard. The average wind field at the area above the boundary layer was strong westerly wind, which can be presumed from faster ground speed bounding east (approximately 370 km h${}^{-\mathrm{1}}$) compared to the ground speed of the west-bound flights (approximately 290 km h${}^{-\mathrm{1}}$). Therefore, it can be inferred that the NO${}_{\mathrm{2}}$ plumes below the aircraft at an altitude of 6000 ft ($\sim$1830 m) tended to propagate northward in alignment with the surface wind field, whereas the plumes at higher altitudes were advected east. This emphasizes the importance of considering the three-dimensional space to elaborate the dispersion and transport pathways of plumes emitted from point sources, as well as accounting for the wind fields at different altitudes. By doing so, the differences between the collocated passive measurements with different viewing geometries, such as satellite, airborne, and ground-based observations, can be better understood.

The NO${}_{\mathrm{2}}$ VCDs over the Jecheon area on 3 November 2020 were slightly lower than the average NO${}_{\mathrm{2}}$ VCDs observed between October and November from 2019 to 2022 (Fig. 6). The NO${}_{\mathrm{2}}$ plumes transported from the SMA (located approximately 100 km to the west) might have affected the area, since the strong westerlies prevailed within both the boundary layer and the free troposphere. However, the effect of the plume was limited since the strong wind must have dispersed the SMA plume in a wide area, lowering the concentration. The uncertainties of the HIS NO${}_{\mathrm{2}}$ VCDs were exceptionally high at Jecheon, with the absolute uncertainty around 0.3 DU, almost 100 % in relative uncertainty. This can be explained by the combined effect of the higher absolute uncertainty level due to low-quality clean pixels used for a SF correction, as well as generally low NO${}_{\mathrm{2}}$ VCDs leading to greater relative uncertainty per same absolute uncertainty. Still, both TROPOMI and HIS were able to capture some increased NO${}_{\mathrm{2}}$ VCDs right above the cement kiln, which has been pointed to as a major industrial point source of NO${}_{\mathrm{2}}$ in the area (Kim et al., 2020).

The research flight over the Pohang area further demonstrated the significance of NO${}_{\mathrm{2}}$ emissions from the steel yards, since one of the world's largest steel yards is located within the city of Pohang (Fig. 7). Southerlies were observed near the surface, and the NO${}_{\mathrm{2}}$ VCDs were higher at the northern side of the steel yard. Since the city of Pohang is much more populated than the cities shown before (i.e., Jecheon, Seosan, Dangjin), the NO${}_{\mathrm{2}}$ VCDs even from the HIS showed peaks over a wider spatial range. Nevertheless, the highest NO${}_{\mathrm{2}}$ VCDs were observed right above the steel yard with a peak concentration of 1.8 DU, which is way beyond the values that can be observed from satellites (i.e., TROPOMI) at the region. For instance, the highest TROPOMI NO${}_{\mathrm{2}}$ VCD was around 0.9 DU. The TROPOMI pixel with the highest VCD was located right next to the exact pixel corresponding to the location of the steel yard and on the northern side. This conforms with southerly dominant atmospheric conditions at the time of observation but is shifted toward the north compared to those found in the airborne HIS observation. The temporal disparities between the HIS frames collected near the steel yard and the TROPOMI overpass were insignificant ($<\mathrm{30}$ min; shown in Fig. S11). Therefore, as an extension of the discussion in Fig. 5g, the conical dispersion and transport of a plume in three-dimensional space must have caused the spatial shift in the peaks due to the different viewing geometries of the satellite and the airborne sensors.

The NO${}_{\mathrm{2}}$ VCDs from the HIS and TROPOMI generally showed good agreement, considering the maximum time difference of 2.5 h between the HIS observation and the TROPOMI overpass, as well as the substantial temporal variability in NO${}_{\mathrm{2}}$ in the vicinity of the emission sources (Park et al., 2022). The correlation coefficient ($r$) ranged from 0.59 to 0.73 (Figs. S9 and S11), except for the research flights at Jecheon with relatively low correlation ($r=\mathrm{0.4}$; Fig. S10). The exceptionally low correlation of the HIS and TROPOMI NO${}_{\mathrm{2}}$ VCDs at Jecheon can be explained by the distinctively high uncertainty of the HIS NO${}_{\mathrm{2}}$ VCDs.

Meanwhile, the HIS NO${}_{\mathrm{2}}$ VCDs exhibited higher levels of uncertainty over the ocean surface (0.10–0.15 DU) than the land surface (0.025–0.075 DU), when the uncertainties associated with the SF calibrations were ruled out. The higher uncertainty observed over the ocean surface can be attributed to the lower surface reflectivity of the ocean (as depicted in Fig. S8). The lower reflectivity can lead to a reduction in the SNR, resulting in an increased level of uncertainty arising from instrumental noise. Furthermore, it can also contribute by itself as an increasing factor of retrieval uncertainty, as discussed in Sect. 3.4 (shown by Figs. 3, S4, and S5).

4.3

Satellite sub-grid variability in NO${}_{\mathrm{2}}$ near the industrial point sources

The NO${}_{\mathrm{2}}$ VCDs from the airborne HIS observations exhibited greater variability than the collocated TROPOMI. In other words, the HIS has shown considerable sub-grid variability within a footprint of TROPOMI. To elaborate on the difference between the HIS and TROPOMI in respect of satellite sub-grid variability, HIS and TROPOMI NO${}_{\mathrm{2}}$ VCDs from all five research flights were collected. Before comparing the collected set of collocated HIS and TROPOMI NO${}_{\mathrm{2}}$ VCDs, bias offsets were incorporated into the HIS NO${}_{\mathrm{2}}$ VCDs per flight (more precisely, per spectral SF). This adjustment aimed to ensure that the mean bias between the HIS and the TROPOMI NO${}_{\mathrm{2}}$ VCDs is zero for a dataset using the same spectral SF because the HIS NO${}_{\mathrm{2}}$ VCDs can show systematic biases depending on which reference NO${}_{\mathrm{2}}$ VCD was used to calculate spectral SF. Therefore, it is necessary to remove the systematic biases incurred from the spectral SF correction to avoid an unintended increase in mean absolute error (MAE) when all the datasets from different flights are combined for comparison.

As a result, Fig. 8a shows a scatter plot comparing the bias-corrected HIS NO${}_{\mathrm{2}}$ VCDs and the collocated TROPOMI data collected from all five research flights conducted for this study. Since the bias-corrected dataset was used for comparison, the mean bias was zero. The overall correlation coefficient ($r$) was 0.73, which was almost equivalent to the highest correlation found in individual flights. The MAE, which can be the indication of the average TROPOMI sub-grid variability incorporating all the research flights, was 0.106 DU. However, considering the random uncertainty of the HIS NO${}_{\mathrm{2}}$ VCD ranging from 0.025 DU (bright land surface) to 0.15 DU (dark ocean surface) under the circumstances that the adequate calibrations were applied (i.e., clean-pixel correction), these random uncertainties posed in the HIS NO${}_{\mathrm{2}}$ VCDs can explain most of the parts of the mean absolute error. Accordingly, it may be inappropriate to argue that the mean absolute error of 0.106 DU represents the overall sub-grid variability at the targeted domains.

Figure 8

(a) Comparison of bias-corrected HIS and TROPOMI NO${}_{\mathrm{2}}$ VCDs with the corresponding HIS NO${}_{\mathrm{2}}$ VCD quantiles (10th, 25th, 75th, 90th) shown over the TROPOMI NO${}_{\mathrm{2}}$ VCD range. (b) Differences between HIS NO${}_{\mathrm{2}}$ VCDs quantile values (maximum–minimum, 51 %–49 %) at corresponding TROPOMI NO${}_{\mathrm{2}}$ VCD range.

[Figure omitted. See PDF]

Figure 8b is a more explicit depiction of the TROPOMI sub-grid variability according to its NO${}_{\mathrm{2}}$ VCDs. From the quantile values of the HIS NO${}_{\mathrm{2}}$ VCDs calculated for every 0.1 DU bin of the TROPOMI NO${}_{\mathrm{2}}$ VCDs, the differences between the quantiles (maximum–minimum, 51st–49th quantiles) according to TROPOMI NO${}_{\mathrm{2}}$ VCD bins are shown. An evident inclination has been identified in the spread of the HIS NO${}_{\mathrm{2}}$ VCD quantile values as the collocated TROPOMI NO${}_{\mathrm{2}}$ VCD increases. Specifically, the difference between 75th and 25th HIS NO${}_{\mathrm{2}}$ VCD quantile values for the TROPOMI footprints with NO${}_{\mathrm{2}}$ VCD below or equal to 0.6 DU yielded 0.15 DU on average, while for the TROPOMI footprints with NO${}_{\mathrm{2}}$ VCDs greater than or equal to 0.8 DU yielded 0.52 DU.

Various artifacts, even besides the retrieval uncertainty, contribute to the variabilities expressed in Fig. 8. Hence, it is challenging to isolate the actual contribution of satellite sub-grid variability. For instance, temporal discrepancies between the HIS frames and the TROPOMI overpass, as well as the combined effect of different viewing geometry and the complex vertical structure of the wind fields, can result in meaningful disparities between the two collocated measurements. Nevertheless, the notably greater HIS NO${}_{\mathrm{2}}$ VCD variability in the TROPOMI footprints with high NO${}_{\mathrm{2}}$ VCDs compared to footprints with relatively low NO${}_{\mathrm{2}}$ VCDs suggests a noticeable increase in satellite sub-grid variability near significant emission sources. Accordingly, the inherent spatial resolution limitations of the satellites can impact their capability to monitor plumes in close proximity to emission sources. To gain a better understanding of the atmospheric conditions near emission hot spots, it is essential to utilize complementary high-resolution observations, such as airborne hyperspectral observations.

5 Summary and conclusions

In this study, we developed the versatile NO${}_{\mathrm{2}}$ VCD retrieval algorithm, so-called the modified wavelength pair (MWP) method, and have applied it to the spectra obtained from the low-cost airborne hyperspectral imaging sensor (HIS) observations at three (Chungnam, Jecheon, and Pohang) industrial areas in South Korea. The MWP method is based on the pseudo-linear relationship between the wavelength pair radiance ratio ($R$) and the NO${}_{\mathrm{2}}$ VCD. The look-up table (LUT) of $R$ values for each wavelength pair was constructed using forward radiative transfer simulations with UVSPEC, which encompassed all possible conditions of the airborne HIS observations. An analytical solution has been derived to calculate NO${}_{\mathrm{2}}$ VCDs using the $R$ values calculated from the airborne-observed spectrum and the corresponding $R$ values in the LUT. This analytical solution allows for the retrieval of three independent NO${}_{\mathrm{2}}$ VCD values from a single HIS spectrum, one for each set of wavelength pairs. The merits of the MWP method are as follows: (1) applicable to the spectra obtained from instruments with relatively low spectral resolution and unstable characteristics, attributed to the lack of precise maintenance or adequate temperature–vibration management; (2) insensitive to broadband spectral tendencies such as from aerosol and surface reflectivity; and (3) computationally competitive from the simplified LUT calculations. For instance, high-resolution chemical transport model simulations are not essentially required.

An analytical derivation of the “true” NO${}_{\mathrm{2}}$ VCD in the MWP method enables feasible error estimation based on sensitivity tests for the LUT. The greatest source of uncertainty was the uncertainty in the spectral scale factor (SF) calibration, which largely depended on the quality of the clean pixel. The uncertainty in the spectral shift calibrations and the instrument noise was the succeeding source of uncertainty, whereas the uncertainties in variables constructing the LUT were relatively minor. The final airborne HIS NO${}_{\mathrm{2}}$ VCD data at a spatial resolution of approximately 400 m $\times$ 400 m (at a flying altitude of 6000 ft, $\sim$1830 m) exhibited uncertainty ranging from 0.025 DU (for bright land surfaces) to 0.15 DU (for dark ocean surfaces), except for the research flight at Jecheon where only a low-quality clean pixel was applicable. However, it should be noted that additional uncertainties could be found in the assumed NO${}_{\mathrm{2}}$ vertical profile, aerosol optical properties, or the effective temperature of NO${}_{\mathrm{2}}$.

The NO${}_{\mathrm{2}}$ plumes emitted from steel yards, both in Chungnam and Pohang, were particularly prominent among the point sources, while stagnant atmospheric conditions with low surface wind speed were a prerequisite for elevated NO${}_{\mathrm{2}}$ concentrations at the proximal area to the emission sources. The overall correlation coefficient ($r$) between the collocated HIS and TROPOMI NO${}_{\mathrm{2}}$ VCDs was 0.73, with a mean absolute error of 0.106 DU. The correlation between the two was limited by temporal disparities between the HIS frames and the TROPOMI overpass, as well as the different observation geometries combined with complex vertical wind fields. There was an apparent increase in the variability in HIS NO${}_{\mathrm{2}}$ VCDs corresponding to the TROPOMI footprints as the TROPOMI NO${}_{\mathrm{2}}$ VCD increased. This implies an amplification of satellite sub-grid variability in NO${}_{\mathrm{2}}$ VCDs near the center of the plume. Therefore, the inherent limitations of satellite spatial resolution can significantly impact its utility in monitoring plumes near emission sources, and it is necessary to employ complementary high-resolution observations. When low-cost, feasible, and compact hyperspectral imagers, such as the HIS used in this study, are utilized to acquire high-spatial-resolution observations, it will lead to an improved understanding of atmospheric conditions near emission point sources through more feasible and frequent observations.

Code and data availability

The airborne HIS NO${}_{\mathrm{2}}$ VCD data are publicly available at 10.7910/DVN/YCZ9JU (Park, 2023). The TROPOMI NO${}_{\mathrm{2}}$ VCD datasets are available at NASA GES DISC (10.5270/S5P-9bnp8q8; ESA and KNMI, 2021). The ERA5 reanalysis datasets are available at 10.24381/cds.bd0915c6 (Hersbach et al., 2023a) and 10.24381/cds.adbb2d47 (Hersbach et al., 2023b). VIIRS BRDF datasets are from 10.5067/VIIRS/VNP43C1.001 (Schaaf et al., 2019). Other datasets such as the raw airborne HIS spectra, the look-up table, ancillary datasets for calibration, and associated codes for their processing may be shared on request to the corresponding author.

The supplement related to this article is available online at: https://doi.org/10.5194/amt-17-197-2024-supplement.

Author contributions

JUP: conceptualization; formal analysis; methodology; visualization; writing – original draft. HJK, JSP, and JC: data curation, funding acquisition, project administration. SSP: conceptualization; methodology; writing – review and editing. KB, JJL, and CKS: data curation, writing – review and editing. SP, KS, and YC: data curation, writing – review and editing. SWK: conceptualization; funding acquisition; writing – review and editing; supervision.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

We would like to express our gratitude to Beom-Keun Seo and Jongho Kim at Hanseo University for their assistance in arranging and providing technical support for the airborne observations.

Financial support

This study was funded by the Fine Particle Research Initiative in East Asia Considering National Differences (FRIEND) through the National Research Foundation of Korea (NRF), funded by the Ministry of Science and ICT (grant no. 2020M3G1A1114615) and NIER research grant (grant no. NIER-SP2022-268). Jin-Soo Park, Hyun-Jae Kim, and Jinsoo Choi were supported by NIER research grant (grant no. NIER–2023-01-01–142). Chang-Keun Song was supported by a grant from the National Institute of Environment Research (NIER), funded by the Ministry of Environment (MOE) of the Republic of Korea (grant no. NIER-2021-03-03-007).

Review statement

This paper was edited by Lok Lamsal and reviewed by two anonymous referees.