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

High-spatial-resolution NO₂ 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₂ emissions are located. An innovative and versatile approach for NO₂ 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₂ plumes emitted from steel yards were particularly conspicuous among the various NO₂ point sources, with peak NO₂ VCDs of 2.0 DU (Dobson unit) at Chungnam and 1.8 DU at Pohang. Typical NO₂ 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₂ VCDs retrieved from the HIS with the MWP method showed a good correlation with the collocated Tropospheric Monitoring Instrument (TROPOMI) data (r=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₂ 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₂ VCD variability (e.g., difference between 25th and 75th quantiles) over the TROPOMI footprints with NO₂ VCD values exceeding 0.8 DU compared to footprints with NO₂ VCD values below 0.6 DU.

1 Introduction

Reactive nitrogen oxides (NO_x), which is a sum of nitrogen dioxide (NO₂) and nitrogen monoxide (NO), are one of the most important trace gases in a polluted atmosphere. 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₃) and aerosol formation (Chan et al., 2010; Seinfeld and Pandis, 2016; Sillman, 1999), both of which are also the major atmospheric pollutants. 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). 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 NO_x must be complemented by remote sensing observations with wider spatial coverage. NO₂ 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₂ from NO_x titration and photolysis during the daytime (Sillman, 1999), observations of NO₂ concentrations can yield the general concentration of NO_x in the atmosphere.

Since the first global NO₂ 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 3.5×5.5 km² at the nadir, while the GOME had a footprint of 40×320 km². However, even up-to-date environmental satellite sensors such as TROPOMI show inherent limitations in spatial resolution to fully capture the NO₂ 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₂ VCDs from hyperspectral imagers has demonstrated its effectiveness in elucidating the spatial variability in NO₂ near emission hot spots (Judd et al., 2019; Tack et al., 2021). However, the highly resolved airborne observations of NO₂ 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₂ 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₂ 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₂ 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=0.84) and TROPOMI (r=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₂ VCD observations can be further highlighted in East Asia, where NO_x emissions are still substantial, and tropospheric O₃ 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₂ 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₂ 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 2¹⁶, with 1392 (spatial columns) ×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 µ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.

Despite these instrument specifications, the HIS shows almost no sensitivity beyond the UV-A region (λ<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₂ are well distinguished, making it suitable for NO₂ VCD retrievals. For instance, Park et al. (2019) retrieved NO₂ VCDs from ground-based observations of the diffuse radiances using the identical instrument (i.e., HIS) and presented its compatibility for NO₂ VCD retrievals by demonstrating a high correlation (r=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₂ VCD retrieval algorithm

Most modern passive remote sensing techniques aimed to observe atmospheric trace gases (e.g., NO₂, O₃, 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₂ 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₃ and NO₂ 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₂ 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₂ VCD retrieval from the airborne HIS spectra using the MWP method.

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₂ absorption exhibits significant high-frequency spectral variability within the VIS band. These pairs are composed of wavelengths that exhibit weak and strong NO₂ 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 λ₁ shows stronger sensitivity to NO₂ concentration, while the longer wavelength λ₂ exhibits relatively weaker NO₂ sensitivity (Fig. 2). In contrast, Type_B pairs show the opposite pattern, with weaker NO₂ sensitivity at the shorter wavelength (λ₁) and stronger sensitivity at the longer wavelength (λ₂). 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, λ₁ and λ₂ (where λ₁<λ₂), is expected to be pseudo-linearly dependent on the NO₂ VCD. Therefore, NO₂ VCDs can be retrieved by comparing the observed R value with the R values simulated from the RTM for a range of possible NO₂ VCDs (and controlling other possible artifacts).

Table 2 Wavelength pairs used in the NO₂ VCD retrieval algorithm (MWP method).

Figure 2 Example of spectral radiance sensitivity to NO₂ 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.

Multiple wavelength pairs are selected within the fitting window to reduce the uncertainty. However, when the NO₂ 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₂ 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₂ 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₂ 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₂ 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₂ 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₂ 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).

$$\begin{array}{}\text{(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.\end{array}$$

Equation (1) shows the relationship between NO₂ 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_A and a_B represent the slope, while b_A and b_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²) 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 λ₂ is a k fold of the value at λ₁, 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 (<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.

$$\begin{array}{}\text{(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.\end{array}$$

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_obs) and the observed radiance ratio of the wavelength pair (R_obs) can also be expressed with the simulated radiance ratio (R_rtm) together with the spectral dependency factor k (Eq. 3).

$$\begin{array}{}\text{(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.\end{array}$$

Since the simulated NO₂ VCDs should be the same regardless of wavelength pair type, the “true” NO₂ VCD (VCD_True) can be expressed as follows (Eq. 4).

$$\begin{array}{}\text{(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}\end{array}$$

From Eq. (4), R_rtm, A can be shown as Eq. (5).

$$\begin{array}{}\text{(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}}}}\end{array}$$

Then the first equation of the simultaneous equations of Eq. (3) can be reformulated by replacing R_rtm, A with Eq. (5), yielding Eq. (6).

$$\begin{array}{}\text{(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}}\end{array}$$

Equations (3) and (6) can be reorganized to Eq. (7) and then be solved to a single equation (Eq. 8).

$$\begin{array}{}\text{(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.\end{array}$$

$$\begin{array}{}\text{(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}}}}\end{array}$$

Equation (8) can be reorganized as an expression of R_rtm, B, yielding Eq. (9).

$$\begin{array}{}\text{(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}\end{array}$$

Then the final equation, Eq. (10), for the “true” unbiased NO₂ VCD is derived only as a function of the observed R values of the Type_A and Type_B wavelength pairs (R_obs, A, R_obs, B) and the coefficients for the first-order polynomials regressed upon the relationship between simulated R values and the NO₂ VCDs (i.e., a_A, b_A, a_B, and b_B in Eq. 1).

$$\begin{array}{}\text{(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}\end{array}$$

As shown in Eq. (10), in which the NO₂ VCD is expressed without the broadband spectral dependency term (k), the NO₂ 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₂ 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 ±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₂ 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₂) for the RTM simulations. The vertical profile of NO₂ can significantly alter the AMF or the effective absorption of NO₂ 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.

As an alternative, NO₂ vertical profiles from the high-resolution (3 km × 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₂ vertical profiles under certain planetary boundary layer height (PBLH) and NO₂ 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 × 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 × 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₂ vertical profiles within each target domain were grouped by the CMAQ-driven PBLH with a ±0.2 km range from the PBLH entries shown in Table 3. The variability (i.e., quantiles) of NO₂ mixing ratios at each altitude level was calculated from vertical profiles grouped under the same PBLH conditions, and the NO₂ vertical profile corresponding to certain a NO₂ VCD was deduced from Levenberg–Marquardt fitting with the quantile value that best fits the targeted column-integrated NO₂ VCD (Fig. S2). It should be noted that the stratospheric NO₂ VCD (i.e., NO₂ 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₂ 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₂ 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₂, such as surface reflectivity and NO₂ 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₂ 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₂ VCD value corresponding to the clean pixel. Whenever a certain NO₂ VCD value is determined for a clean pixel, it directly affects all the other NO₂ 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₂ VCD, which may be biased proportionally to the bias of NO₂ 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₂ VCDs retrieved from the MWP method can be estimated from the errors in input parameters and their sensitivity to the resulting NO₂ VCD. One of the advantages of the MWP method is that the solution (NO₂ VCD) is retrieved analytically, making the uncertainty estimation feasible by adopting Gaussian error propagation. The error estimation of the NO₂ 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_A, b_A, a_B, and b_B in Eq. 10) are calculated from the NO₂ 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₂ 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₂ VCD retrieval in the MWP method.

To calculate error propagation, R_obs values in Eq. (10) should be transposed with R_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₂ VCD. Furthermore, we can transpose the ratio between R_rtm, A and R_rtm, 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_rtm, A and R_rtm, B for the error propagation.

$$\begin{array}{}\text{(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}\end{array}$$

Under these circumstances, the error in VCD_True can be estimated from Eq. (11) as follows (Eq. 12).

$$\begin{array}{}\text{(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}}\end{array}$$

The uncertainty of the Q value (σ_Q) is determined by the sensitivity of Q values to the uncertainty of input parameters. However, σ_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 σ_Q by multiplying the Q value calculated from the HIS spectrum in the actual NO₂ VCD retrieval.

$$\begin{array}{}\text{(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)\end{array}$$

Table 4 presents the assumed uncertainty range of input parameters (σ_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 10⁻² and 10⁻¹ 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 (σ_var) assumed for the error estimation.

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₂ 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).

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 2×10³, 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₂ VCD conditions. NO₂ 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₂ 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₂ 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₂ vertical profiles, as well as consequential structural uncertainty caused by varying effective temperature of NO₂ (i.e., temperature dependency of NO₂ 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₂ 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₂ 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₂ 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 <0.3) TROPOMI NO₂ 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₂ 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₂ VCDs were observed in all three target domains. The NO₂ 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₂ 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₂ VCD_Total in October and November from 2019 to 2022. TROPOMI swath datasets were downscaled to grids of 0.25^∘ × 0.25^∘ for South Korea and 0.01^∘ × 0.01^∘ for the target domains (i.e., Chungnam, Jecheon, and Pohang) by calculating average NO₂ 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.

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 (∼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^∘), resulted in a narrow observation swath ranging between 340 m (at 5000 ft, ∼1520 m) and 750 m (at 11 000 ft, ∼3350 m).

4.2 NO₂ VCDs in three industrial areas

The MWP method applied to all three sets of wavelength pairs yields three independent NO₂ VCDs with corresponding uncertainty estimates from a single HIS spectrum. To take full advantage of three independent NO₂ VCDs, the weighted mean value was taken as a single representative NO₂ 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₂ VCD from three different wavelength pair sets (Eq. 14), which can minimize the random uncertainty theoretically (Xie et al., 2001).

$$\begin{array}{}\text{(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)}}\end{array}$$

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, ∼3350 m, and 13 m at 5000 ft, ∼1520 m) compared to the along-track length scale (approximately 100 m under an aircraft ground speed of 360 km h⁻¹ and 1 s of exposure). Moreover, the remaining uncertainties of the NO₂ 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₂ 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₂ VCD (as shown in Eq. 14). Therefore, the final HIS data have a footprint almost symmetric at an altitude of 6000 ft (∼1830 m; 400 m × 400 m) but can be variable according to the ground speed of an airplane and the flying altitude.

The NO₂ 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₂ VCDs observed from three research flights over the Chungnam area showed different spatial patterns due to different wind fields (Fig. 5). The highest NO₂ concentration was observed on 24 November 2022 when the wind was calm (<2 m s⁻¹) within the boundary layer (Fig. 5e). Since the atmosphere was stagnant near the surface, the peaks of NO₂ VCDs were observed right above the major emission point sources. The highest concentration was found over the steel yard, with the peak HIS NO₂ 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₂ VCDs over the petrochemical complex (>1.2 DU), while the NO₂ VCD of the corresponding TROPOMI footprint remained relatively low (0.9 DU). While the small elevated NO₂ 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₂ 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₂ and populated urban areas denoted with red and blue circles, respectively. Panels (c), (e), and (g) show NO₂ VCD retrieval results from HIS overlaid on TROPOMI swath data, and panels (d), (f), and (h) show the uncertainty of HIS-retrieved NO₂ 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 6 (a) Downscaled average TROPOMI NO₂ 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₂ 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₂ VCD retrieval results from HIS overlaid on TROPOMI swath data, and panel (d) shows the uncertainty of HIS-retrieved NO₂ VCD.

Figure 7 (a) Downscaled average TROPOMI NO₂ 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₂ 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₂ VCD retrieval results from HIS overlaid on TROPOMI swath data, and panel (d) shows the uncertainty of HIS-retrieved NO₂ VCD.

The research flights on 17 October 2020 (Fig. 5c) and 25 November 2022 (Fig. 5g) exhibited relatively lower NO₂ 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₂ plume conically dispersed and transported to the north under clear southerly wind at the boundary layer. The significance of NO₂ 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₂ 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₂ 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₂ 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⁻¹) compared to the ground speed of the west-bound flights (approximately 290 km h⁻¹). Therefore, it can be inferred that the NO₂ plumes below the aircraft at an altitude of 6000 ft (∼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₂ VCDs over the Jecheon area on 3 November 2020 were slightly lower than the average NO₂ VCDs observed between October and November from 2019 to 2022 (Fig. 6). The NO₂ 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₂ 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₂ VCDs leading to greater relative uncertainty per same absolute uncertainty. Still, both TROPOMI and HIS were able to capture some increased NO₂ VCDs right above the cement kiln, which has been pointed to as a major industrial point source of NO₂ in the area (Kim et al., 2020).

The research flight over the Pohang area further demonstrated the significance of NO₂ 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₂ 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₂ VCDs even from the HIS showed peaks over a wider spatial range. Nevertheless, the highest NO₂ 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₂ 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 (<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₂ 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₂ 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=0.4; Fig. S10). The exceptionally low correlation of the HIS and TROPOMI NO₂ VCDs at Jecheon can be explained by the distinctively high uncertainty of the HIS NO₂ VCDs.

Meanwhile, the HIS NO₂ 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₂ near the industrial point sources

The NO₂ 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₂ VCDs from all five research flights were collected. Before comparing the collected set of collocated HIS and TROPOMI NO₂ VCDs, bias offsets were incorporated into the HIS NO₂ VCDs per flight (more precisely, per spectral SF). This adjustment aimed to ensure that the mean bias between the HIS and the TROPOMI NO₂ VCDs is zero for a dataset using the same spectral SF because the HIS NO₂ VCDs can show systematic biases depending on which reference NO₂ 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₂ 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₂ 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₂ 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₂ VCDs with the corresponding HIS NO₂ VCD quantiles (10th, 25th, 75th, 90th) shown over the TROPOMI NO₂ VCD range. (b) Differences between HIS NO₂ VCDs quantile values (maximum–minimum, 51 %–49 %) at corresponding TROPOMI NO₂ VCD range.

Figure 8b is a more explicit depiction of the TROPOMI sub-grid variability according to its NO₂ VCDs. From the quantile values of the HIS NO₂ VCDs calculated for every 0.1 DU bin of the TROPOMI NO₂ VCDs, the differences between the quantiles (maximum–minimum, 51st–49th quantiles) according to TROPOMI NO₂ VCD bins are shown. An evident inclination has been identified in the spread of the HIS NO₂ VCD quantile values as the collocated TROPOMI NO₂ VCD increases. Specifically, the difference between 75th and 25th HIS NO₂ VCD quantile values for the TROPOMI footprints with NO₂ VCD below or equal to 0.6 DU yielded 0.15 DU on average, while for the TROPOMI footprints with NO₂ 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₂ VCD variability in the TROPOMI footprints with high NO₂ VCDs compared to footprints with relatively low NO₂ 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₂ 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₂ 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₂ 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₂ 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₂ 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₂ VCD data at a spatial resolution of approximately 400 m × 400 m (at a flying altitude of 6000 ft, ∼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₂ vertical profile, aerosol optical properties, or the effective temperature of NO₂.

The NO₂ 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₂ concentrations at the proximal area to the emission sources. The overall correlation coefficient (r) between the collocated HIS and TROPOMI NO₂ 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₂ VCDs corresponding to the TROPOMI footprints as the TROPOMI NO₂ VCD increased. This implies an amplification of satellite sub-grid variability in NO₂ 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.