Tropospheric NO₂ measurements using a three-wavelength optical parametric oscillator differential absorption lidar

The conventional two-wavelength differential absorption lidar (DIAL) has measured air pollutants such as nitrogen dioxide (NO₂). However, high concentrations of aerosol within the planetary boundary layer (PBL) can cause significant retrieval errors using only a two-wavelength DIAL technique to measure NO₂. We proposed a new technique to obtain more accurate measurements of NO₂ using a three-wavelength DIAL technique based on an optical parametric oscillator (OPO) laser. This study derives the three-wavelength DIAL retrieval equations necessary to retrieve vertical profiles of NO₂ in the troposphere. Additionally, two rules to obtain the optimum choice of the three wavelengths applied in the retrieval are designed to help increase the differences in the NO₂ absorption cross-sections and reduce aerosol interference. NO₂ retrieval relative uncertainties caused by aerosol extinction, molecular extinction, absorption of gases other than the gas of interest and backscattering are calculated using two-wavelength DIAL (438 and 439.5 nm) and three-wavelength DIAL (438, 439.5 and 441 nm) techniques. The retrieval uncertainties in aerosol extinction using the three-wavelength DIAL technique are reduced to less than 2 % of those when using the two-wavelength DIAL technique. Moreover, the retrieval uncertainty analysis indicates that the three-wavelength DIAL technique can reduce more fluctuation caused by aerosol backscattering than the two-wavelength DIAL technique. This study presents NO₂ concentration profiles which were obtained using the HU (Hampton University) three-wavelength OPO DIAL. As a first step to assess the accuracy of the HU lidar NO₂ profiles, we compared the NO₂ profiles to simulated data from the Weather Research and Forecasting Chemistry (WRF-Chem) model. This comparison suggests that the NO₂ profiles retrieved with the three-wavelength DIAL technique have similar vertical structure and magnitudes typically within ±0.1 ppb compared to modeled profiles.

1 Introduction

Nitrogen dioxide (NO₂) plays a critical role in the tropospheric chemistry and is one of the reactive gases collectively referred to as “nitrogen oxides” (NO_x = nitric oxide and nitrogen dioxide, NO + NO₂) (U.S. EPA, 2020). The sources of NO_x emissions include transportation (on-road vehicles, airplanes, trains, ships), wood burning, industrial and chemical processes, activities for oil and gas development, soil emissions, lightning, and wildfires (see Nitrogen Oxides Emissions indicator) (U.S. EPA, 2020). Once emitted, NO reacts rapidly in the presence of ozone to form NO₂. In US urban locations, most measured airborne NO₂ comes from the reaction of these two precursors rather than from direct NO₂ emissions (Bertram, et al., 2005; Beirle, et al., 2011). Scientific evidence indicates that short-term NO₂ exposure, ranging from 30 min to 24 h, can cause the exacerbation of asthma symptoms, in some cases resulting in hospitalization (Berglund, et al., 1993). Long-term NO₂ exposure is likely to have a causal relationship with respiratory effects, based on evidence for the development of asthma (U.S. EPA, 2016). And NO₂ will be included in future cycles of the Global Burden of Disease as global exposure estimates and evidence of their role as independent risk factors accumulate (Larkin et al., 2017). Additionally, atmospheric processing of NO₂ leads to the formation of nitrogen-bearing particles that can eventually deposit to the surface, causing acidification, nitrogen enrichment and other ecological effects (Russell et al., 2012). Local or global NO₂ monitoring is essential for understanding atmospheric chemistry as well as for human-health and environmental management and control.

Measurements of the intensity of ultraviolet or visible absorption spectra from the ground or from satellites are commonly used to retrieve the column density of NO₂ (Celarier et al., 2008; Valks et al., 2011; Berg et al., 2012). Satellite-based instruments such as the Ozone Monitoring Instrument (OMI), Global Ozone Monitoring Experiment (GOME and GOME-2) and SCanning Imaging Absorption SpectroMeter for Atmospheric CHartographY (SCIAMACHY) can provide global-scale NO₂ column measurements during daytime (Boersma et al., 2008; Bucsela et al., 2008). Moreover, plumes of NO₂ by cities, power plants and even ships have been able to be tracked using the recent high-spatial-resolution observations of NO₂ from the TROPOspheric Monitoring Instrument (TROPOMI) on Sentinel-5P since 2017 (Lorente, et al., 2019; Georgoulias et al., 2020). However, they are unable to obtain local high-temporal-resolution NO₂ emissions such as variations in hourly NO₂ concentrations due to their long repeat cycle since the lifetime of tropospheric NO₂ is only about 6 h in summer and 18–24 h in winter due to the photochemical effect (Beirle, et al., 2003; Cui et al., 2016). In addition, measurements of tropospheric NO₂ from satellites or aircraft are also influenced and limited by clouds (Bovensmann et al., 1999; Liang et al., 2017). Ground-based measurements of column NO₂ from instruments such as Pandora using differential optical absorption spectroscopy (DOAS) are often used for the validation of satellite instruments (Herman et al., 2009; Lamsal et al., 2014; Kollonige et al., 2018). In situ measurements of near-surface NO₂ can best monitor local emissions. However, at this point in time, they cannot provide vertically resolved measurements. Balloon measurements using a NO₂ sonde can produce vertical profiles, but these measurements are very limited in time and space, especially in the Southern Hemisphere. The primary source of data on the vertical distribution of NO₂ comes from operational sites around the world. However, their operation can be expensive and labor-intensive (Scott et al., 1999; Herman et al., 2009; Sluis et al., 2010).

The differential absorption lidar (DIAL) technique offers the potential for autonomous, 24 × 7 operation, with improved temporal resolution. Absorption of light by molecules is the basis for DIAL and numerous atmospheric constituents absorbing light. Conventional DIAL operates at two absorbing wavelengths, one stronger than the other, indicated by the on (λ_on) and off (λ_off) wavelength of the gaseous absorption feature of interest. Because of different absorption at λ_on and λ_off, the difference between the backscattered laser signals at the two wavelengths can be used to derive the number density of the absorption gas. Taking the log-ratio of these returns at closely spaced wavelengths removes system parameters and attenuation to and from the target of interest (Rothe et al., 1974; Sullivan et al., 2014). Thus, this technology provides measurements of the concentration of gas, such as NO₂, O₃ and SO₂, at a particular location and time (Fredriksson and Hertz, 1984; Newchurch et al., 2003; Volten et al., 2009; Kuang et al., 2013; Sullivan et al., 2017). The DIAL technique provides the unique capability of remotely monitoring localized NO₂ concentrations and emissions in urban and rural areas and profiling their tropospheric vertical NO₂ concentration. However, aerosols are abundant within the planetary boundary layer (PBL) and can cause significant retrieval errors in a two-wavelength DIAL technique to measure NO₂. To better understand this aerosol problem and produce a more accurate NO₂ profile measurement, we described a new technique using a three-wavelength DIAL technique based on the intrinsic capabilities of using a multi-wavelength optical parametric oscillator (OPO) laser system. Hampton University (HU) has incorporated an OPO laser into its lidar system. The OPO laser enables researchers to optimize (tune) wavelength choices for specific measurements (Weibring et al., 2003). The three-wavelength DIAL retrieval equations are derived in this study. Our optimum choices for the three wavelengths to be used for our NO₂ retrievals are designed to help increase the difference in NO₂ absorption cross-section and reduce aerosol influence. NO₂ retrieval relative uncertainties are calculated using the two-wavelength DIAL (438 and 439.5 nm) and the three-wavelength DIAL (438, 439.5 and 441 nm). Tropospheric NO₂ profiles were obtained by applying the proposed technique to HU OPO DIAL lidar. As a first-order assessment, the HU lidar results were compared with simulated data from the Weather Research and Forecasting Chemistry (WRF-Chem) air quality model.

2 Method

To minimize aerosols-interference in the retrievals of NO₂, a three-wavelength DIAL technique was proposed with ${\mathit{\lambda }}_{\mathrm{1}}<{\mathit{\lambda }}_{\mathrm{2}}<{\mathit{\lambda }}_{\mathrm{3}}$. Table 1 shows expressions for the extinction and backscatter of molecules and aerosols for these three wavelengths. In Table 1, β_m and β_a are backscatter from molecules and aerosols for the wavelength of λ₂; α_m and α_a are the extinction of molecules and aerosols for the wavelength of λ₂; and e is the aerosol Ångström exponent, assumed to be equal for the three wavelengths because the three wavelengths are very close.

Table 1 Extinction and backscatter of molecule and aerosol for wavelengths of λ₁, λ₂ and λ₃.

The three elastic lidar equations can be expressed as

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}& X\left({\mathit{\lambda }}_{\mathrm{1}},Z\right)=C_{\mathrm{1}}{\displaystyle \frac{\left[{\left(\frac{{\mathit{\lambda }}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{2}}}\right)}^{-\mathrm{4}}{\mathit{\beta }}_{\mathrm{m}}\left(Z\right)+{\left(\frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{1}}}\right)}^{-e}{\mathit{\beta }}_{\mathrm{a}}\left(Z\right)\right]}{Z^{\mathrm{2}}}}\\ & \times \exp \mathit{\{}-\mathrm{2}\underset{\mathrm{0}}{\overset{Z}{\int }}[{\left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{2}}}}\right)}^{-\mathrm{4}}{\mathit{\alpha }}_{\mathrm{m}}\left(Z\right)+{\left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{1}}}}\right)}^{-e}{\mathit{\alpha }}_{\mathrm{a}}\left(Z\right)\\ & +{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{1}},Z\right)N_N\left(Z\right)+O_{\mathrm{abs}}\left({\mathit{\lambda }}_{\mathrm{1}},Z\right)\left]\mathrm{d}z\mathit{\right\}}\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}& X\left({\mathit{\lambda }}_{\mathrm{2}},Z\right)=C_{\mathrm{2}}{\displaystyle \frac{\left[{\mathit{\beta }}_{\mathrm{m}}\left(Z\right)+{\mathit{\beta }}_{\mathrm{a}}\left(Z\right)\right]}{Z^{\mathrm{2}}}}\exp \mathit{\{}-\mathrm{2}\underset{\mathrm{0}}{\overset{Z}{\int }}\left[{\mathit{\alpha }}_{\mathrm{m}}\left(z\right)\right.\\ & \left.+{\mathit{\alpha }}_{\mathrm{a}}\left(z\right)+{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{2}},z\right)N_N\left(z\right)+O_{\mathrm{abs}}\left({\mathit{\lambda }}_{\mathrm{2}},z\right)\right]\mathrm{d}z\mathit{\}}\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}& X\left({\mathit{\lambda }}_{\mathrm{3}},Z\right)=C_{\mathrm{3}}{\displaystyle \frac{\left[{\left(\frac{{\mathit{\lambda }}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{2}}}\right)}^{-\mathrm{4}}{\mathit{\beta }}_{\mathrm{m}}\left(Z\right)+{\left(\frac{{\mathit{\lambda }}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{2}}}\right)}^{-e}{\mathit{\beta }}_{\mathrm{a}}\left(Z\right)\right]}{Z^{\mathrm{2}}}}\\ & \times \exp \mathit{\{}-\mathrm{2}\underset{\mathrm{0}}{\overset{Z}{\int }}[{\left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{2}}}}\right)}^{-\mathrm{4}}{\mathit{\alpha }}_{\mathrm{m}}\left(z\right)+{\left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{2}}}}\right)}^{-e}{\mathit{\alpha }}_{\mathrm{a}}\left(z\right)\\ & +{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{3}},z\right)N_N\left(z\right)+O_{\mathrm{abs}}\left({\mathit{\lambda }}_{\mathrm{3}},z\right)\mathrm{d}z\left]\mathit{\right\}},\end{array}\end{array}$$

where X is the lidar signal; C₁, C₂ and C₃ are lidar constants; the subscripts a and m represent aerosol and molecule, respectively; σ_N is the absorption cross-section for the gas of interest; N_N is the molecular density of the gas of interest; O_abs is absorption of gases other than the gas of interest; and z is the altitude. The molecular density of the gas of interest can be obtained using Eqs. (1), (2) and (3).

The NO₂ density retrieval equation can be expressed as

$$\begin{array}{}\text{(4)}& {\displaystyle }{N}_N\left(Z\right)={\displaystyle \frac{\begin{array}{c}\frac{\mathrm{1}}{\mathrm{2}}\times \frac{\mathrm{d}}{\mathrm{d}z}\left[\ln \frac{X\left({\mathit{\lambda }}_{\mathrm{1}},Z\right)X\left({\mathit{\lambda }}_{\mathrm{3}},Z\right)}{X{\left({\mathit{\lambda }}_{\mathrm{2}},Z\right)}^{\mathrm{2}}}\right]-\mathrm{AED}\left(z\right)\\ -\mathrm{MED}\left(z\right)-\mathrm{OAD}\left(z\right)-B\left(z\right)\end{array}}{\mathrm{\Delta }{\mathit{\sigma }}_N}}\text{(5)}& {\displaystyle }\mathrm{\Delta }{\mathit{\sigma }}_N=\mathrm{2}{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{2}}\right)-{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{1}}\right)-{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{3}}\right)\text{(6)}& {\displaystyle }B\left(z\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left[\ln {\displaystyle \frac{\begin{array}{c}\left[{\left(\frac{{\mathit{\lambda }}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{2}}}\right)}^{-\mathrm{4}}{\mathit{\beta }}_{\mathrm{m}}\left(Z\right)+{\left(\frac{{\mathit{\lambda }}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{2}}}\right)}^{-e}{\mathit{\beta }}_{\mathrm{a}}\left(Z\right)\right]\\ \times \left[{\left(\frac{{\mathit{\lambda }}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{2}}}\right)}^{-\mathrm{4}}{\mathit{\beta }}_{\mathrm{m}}\left(Z\right)+{\left(\frac{{\mathit{\lambda }}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{2}}}\right)}^{-e}{\mathit{\beta }}_{\mathrm{a}}\left(Z\right)\right]\end{array}}{{\left[{\mathit{\beta }}_{\mathrm{m}}\left(Z\right)+{\mathit{\beta }}_{\mathrm{a}}\left(Z\right)\right]}^{\mathrm{2}}}}\right]\text{(7)}& {\displaystyle }\mathrm{AED}\left(z\right)=K{\mathit{\alpha }}_{\mathrm{a}}\left(Z\right)K=\mathrm{2}-{\left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{2}}}}\right)}^{-e}-{\left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{2}}}}\right)}^{-e}\text{(8)}& {\displaystyle }\mathrm{MED}\left(z\right)=\left[\mathrm{2}-{\left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{2}}}}\right)}^{-\mathrm{4}}-{\left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{2}}}}\right)}^{-\mathrm{4}}\right]{\mathit{\alpha }}_{\mathrm{m}}\left(Z\right)\text{(9)}& {\displaystyle }\mathrm{OAD}\left(z\right)=\mathrm{2}{O}_{\mathrm{abs}}\left({\mathit{\lambda }}_{\mathrm{2}},z\right)-O_{\mathrm{abs}}\left({\mathit{\lambda }}_{\mathrm{1}},z\right)-O_{\mathrm{abs}}\left({\mathit{\lambda }}_{\mathrm{3}},z\right),\end{array}$$

where AED, MED, OAD and B are the correction terms of aerosol extinction, molecular extinction, absorption of gases other than the gas of interest and backscattering, respectively. Because the atmospheric molecular density is relatively stable, MED can be corrected using a numerical model or local real-time radiosonde data. OAD can be removed by choosing appropriate wavelengths. However, aerosol is variable, especially in the PBL. For correction of AED and B, we need accurate aerosol measurements. However, accurate aerosol measurements are not easily obtained. From the above NO₂ retrieval relative equation, all of the correction terms are related to the three wavelengths, so how to choose the three wavelengths is very critical to reduce correction terms and improve the accuracy of NO₂ retrievals. We designed two rules to obtain the optimum choice for the three wavelengths:

• a.

  The three chosen wavelengths increase differences in the NO₂ absorption cross-section (Δσ_N) to improve NO₂ retrieval.

  According to Eq. (4), the greater Δσ_N is, the less all of the correction terms are. So the three chosen wavelengths should help to increase Δσ_N. Generally, researchers only used an increasing absorption method (${\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{1}}\right)<{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{2}}\right)<{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{3}}\right)$) or a decreasing absorption method (${\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{1}}\right)>{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{2}}\right)>{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{3}}\right)$) (illustrated in Fig. 1) to choose the three wavelengths (Wang, et al., 1997; Liu, et al., 2017). Wang used three wavelengths corresponding to the strong, medium and weak absorption of O₃ to obtain an accurate stratospheric ozone profile in the presence of volcanic aerosols. Liu used three wavelengths of 448.10, 447.20 and 446.60 nm corresponding to the strong, medium and weak absorption of NO₂ to retrieve NO₂. Equations (10) and (11) are calculated values of Δσ_N for the increasing absorption method and the decreasing absorption method using Eq. (5). Using the increasing absorption method and the decreasing absorption method to choose the three wavelengths, the values of Δσ_N are both decreased according to Eqs. (10) and (11) compared to the conventional two-wavelength DIAL technique (the conventional two-wavelength DIAL technique retrieval equation is in the Supplement). According to characteristics of the NO₂ absorption spectrum shown in Fig. 2, a bumping absorption method (σ_N(λ₁)<σ_N(λ₂) & σ_N(λ₃)<σ_N(λ₂)) is designed to choose the three wavelengths which can increase the value of Δσ_N compared to the two-wavelength DIAL technique according to Eq. (12). However, for DIAL systems to measure other atmospheric gases like ozone, it is only practical to use wavelength selection Method B because of the shape of the ozone absorption spectrum (lacking narrow peaks).

  Increasing absorption method:

  $$\begin{array}{}\text{(10)}& \begin{array}{rl}\mathrm{\Delta }{\mathit{\sigma }}_N& =\mathrm{abs}\left[{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{2}}\right)-{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{1}}\right)\right]\\ & -\mathrm{abs}\left[{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{2}}\right)-{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{3}}\right)\right]\end{array}\end{array}$$

  Decreasing absorption method:

  $$\begin{array}{}\text{(11)}& \begin{array}{rl}\mathrm{\Delta }{\mathit{\sigma }}_N& =\mathrm{abs}\left[{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{2}}\right)-{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{3}}\right)\right]\\ & -\mathrm{abs}\left[{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{2}}\right)-{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{1}}\right)\right]\end{array}\end{array}$$

  Bumping absorption method:

  $$\begin{array}{}\text{(12)}& \begin{array}{rl}\mathrm{\Delta }{\mathit{\sigma }}_N& =\mathrm{abs}\left[{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{2}}\right)-{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{1}}\right)\right]\\ & +\mathrm{abs}\left[{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{2}}\right)-{\mathit{\sigma }}_N\left({\mathit{\lambda }}_{\mathrm{3}}\right)\right]\end{array}\end{array}$$
• b.

  The three chosen wavelengths can reduce or remove AED.

  It means the value of AED is equal or close to 0. Choosing the appropriate three wavelengths to make the value of K in Eq. (12) equal or close to 0, the value of AED will be equal or close to 0. The value of K in Eq. (12) changes with different aerosol Ångström exponents. For example, to remove boundary layer aerosol influence, we can set aerosol Ångström exponents to 1 to calculate the value of K to choose the three wavelengths because the size of aerosol in the boundary layer is typically large (Schuster et al., 2006; Li et al., 2017).

Figure 1 The chosen three-wavelength methods: increasing absorption method (a), decreasing absorption method (b) and bumping absorption method (c).

Figure 2 NO₂ strong-absorption cross-section between 420 and 450 nm.

Figure 3 Continuum Horizon II energy outputs (a) and parameters (b) with PL 8000 pump.

3 HU three-wavelength OPO DIAL system

The HU lidar is located on the campus of HU (37.02^∘ N, 76.34^∘ W) in Hampton, VA. A Continuum Horizon II tunable OPO laser and a Continuum Powerlite DLS 8000 pump laser have recently been incorporated into the HU lidar system. The OPO laser enables researchers to optimize (tune) the wavelength choices and provides more flexibility than fixed-frequency wavelength shifters such as Raman cells. The wavelength tuning range of our OPO extends from 192 to 2750 nm. This range is fully automated with precision scanning for true hands-free operation. Figure 3a and b show the Continuum Horizon II output energy and its parameters. The OPO laser energy outputs between 400 and 500 nm, which overlaps with the NO₂ strong-absorption spectral zone in Fig. 2 produced near the maximum possible power in the spectrum. Combining the OPO laser energy outputs, NO₂ absorption spectrum and two chosen three-wavelength rules, 438, 439.5 and 441 nm shown in Fig. 2 result in the wavelengths of the HU three-wavelength DIAL system because Δσ_N of the three-wavelength pair is more than other three-wavelength pairs in NO₂ strong-absorption spectral zone, and the K value of the three-wavelength DIAL system is 0.000023 (close to 0). The HU lidar system currently consists of a Continuum OPO laser system as the light source, a 48 in. (1.22 m) non-coaxial Cassegrainian-configured telescope receiver, a light separation system that uses beam splitters and interference filters, a detecting system including photomultiplier tubes (PMTs) and avalanche photodiodes (APDs), and a Licel optical transient recorder. A schematic of the lidar system is shown in Fig. 4. The system can be configured to measure multi-wavelength aerosols and NO₂ density. High-resolution backscatter measurements extend from the boundary layer (1.2 km) to the free troposphere. The pump laser operates at three fixed wavelengths (1064, 532 and 354.7 nm). The 354.7 nm laser is mostly reflected into the OPO laser to produce three wavelengths (338, 339.5 and 441 nm). Steering mirrors whose axes are aligned with a receiving telescope axis direct these laser outputs into the atmosphere. The laser backscatter is collected by a 48 in. (1.22 m) diameter telescope and split into specific wavelength bands by a beam separation unit, which combines filters and beam splitters for dispersion of the return backscatter to various detection channels. Using filters and beam splitters makes the beam-splitting system simple, compact, and easy to change or add other spectral channels for other measurements. Currently, wavelengths of 438, 439.5, 441, 354.7, 532 and 1064 nm are focused on PMTs and APDs and recorded by a Licel data-collecting system for measurements of aerosol and NO₂.

Figure 4 HU lidar system (L: lens; M: mirror; BS: beam splitter; IF: interference filter; FOV: field of view; PMT: photomultiplier tube; APD: avalanche photodetector).

Figure 5 Atmospheric profiles used for modeling NO₂ lidar correction terms. (a) Aerosol extinction profile (red) at 532 nm measured by the HU lidar and molecular extinction profile (blue) at 532 nm derived from local radiosonde data, (b) NO₂ (red) and O₃ (blue) mixing ratio profiles from NASA DISCOVER-AQ and TOLNet.

To demonstrate that the HU three-wavelength OPO DIAL system can effectively reduce aerosol influence and accurately retrieve NO₂, retrieval correction terms of AED, MED, OAD and B in Eq. (4) are simulated using the two-wavelength DIAL technique (438 and 439.5 nm) and the three-wavelength DIAL technique (438, 439.5 and 441 nm). Ozone was used for the simulation of OAD because only ozone absorption can produce a little influence on NO₂ retrieval based on the HITRAN 1.1.2.0 database. Atmospheric data of aerosol, molecules, O₃ and NO₂ for these simulations are from the HU local lidar aerosol measurements, radiosonde, NASA Tropospheric Ozone Lidar Network (TOLNet), and NASA Deriving Information on Surface Conditions from COlumn and VERtically Resolved Observations Relevant to Air Quality (DISCOVER-AQ) measurements shown in Fig. 5. Extinction and backscatter of aerosol at 438, 439.5 and 441 nm can be calculated from the aerosol extinction profile at 532 nm in Fig. 5a with the setting of the lidar ratio = 50 and e=1, 2 and 3. The lidar ratio is wavelength-dependent, and its value in the visible band is in general smaller than in the UV band for the same type of aerosols (Kuang et al., 2020; Reid et al., 2017). Absorption of NO₂ and O₃ at 438, 439.5 and 441 nm can be calculated using their mixing ratio profiles in Fig. 5b and their absorption cross-sections from the HITRAN 1.1.2.0 database. MED, AED, OAD, B and absorption difference in NO₂ (NAD) are simulated using the two-wavelength DIAL technique with different aerosol Ångström exponents (e=1, 2 and 3) shown in Fig. 6a, c and e and the three-wavelength DIAL technique shown in Fig. 6b, d and f. In Fig. 6, red lines are NAD, black lines are MED, deep-blue lines are AED, and light-blue lines are OAD. In Fig. 6, all OAD is far less than NAD. It is concluded that ozone absorption has a negligible influence on the retrieval of NO₂. In Fig. 6a, c and d, MED and AED in the PBL are both more than NAD using the two-wavelength DIAL technique. Because atmospheric molecules are relatively stable, MED can be corrected using local model or real-time radiosonde data. However, aerosol is variable, so aerosols are a significant uncertainty for retrieving NO₂ with the conventional two-wavelength DIAL technique. In Fig. 6b, d and f, MED and AED in the boundary layer are both much smaller than NAD using the proposed three-wavelength DIAL technique. It is proven that the three-wavelength DIAL technique can effectively decrease retrieval errors caused by aerosol extinction. From Fig. 5, we can see that AED using the three-wavelength DIAL technique can be reduced to less than 2 % of AED using the two-wavelength DIAL technique at least. Therefore, even if AED is not corrected, NO₂ can still be accurately retrieved. Moreover, simulated B using the two-wavelength DIAL technique and the three-wavelength DIAL technique is shown in Fig. 6 with green lines. The sharp change in vertical adjacent aerosol backscatter can cause drastic changes in the B term. In Fig. 6, the value of the B term using the three-wavelength DIAL technique is far less than using the two-wavelength DIAL technique. So the three-wavelength DIAL technique can reduce more fluctuation caused by aerosol backscattering than the two-wavelength DIAL technique.

Figure 6 Simulated MED (black), NAD (red), OAD (light blue), AED (deep blue) and B (green) using the two-wavelength DIAL technique with e = 1 (a), e = 2 (c) and e = 3 (e) and the three-wavelength DIAL technique with e = 1 (b), e = 2 (d) and e = 3 (f).

Figure 7 NO₂ number density relative uncertainty owing to air number density (a) and ozone number density (b).

Figure 8 NO₂ number density relative uncertainty owing to aerosol extinction (a) and backscatter (b) with e = 1.

Figure 9 NO₂ number density relative uncertainty owing to aerosol extinction (a) and backscatter (b) with e = 2.

Figure 10 NO₂ number density relative uncertainty owing to aerosol extinction (a) and backscatter (b) with e = 3.

4 Uncertainty analysis

According to NO₂ retrieval Eq. (4), the NO₂ measurement uncertainty is due to molecules, absorption of gases other than the gas of interest, aerosol and noise of lidar signals. The total relative uncertainty can be expressed as Eq. (13) (Leblanc et al., 2016).

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& U_{{\mathrm{NO}}_{\mathrm{2}}}\left(z\right)=\\ & \sqrt{U_{\mathrm{AED}}(z{)}^{\mathrm{2}}+U_{\mathrm{MED}}(z{)}^{\mathrm{2}}+U_{\mathrm{OAD}}(z{)}^{\mathrm{2}}+U_B(z{)}^{\mathrm{2}}+U_{\mathrm{S}}(z{)}^{\mathrm{2}}}\end{array}\end{array}$$

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{U}_{\mathrm{MED}}\left(z\right)& ={\displaystyle \frac{u\left[\mathrm{MED}\left(z\right)\right]}{N_N\left(z\right)\mathrm{\Delta }{\mathit{\sigma }}_N}}\\ & ={\displaystyle \frac{\left[\mathrm{2}-{\left(\frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{1}}}\right)}^{\mathrm{4}}-{\left(\frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{3}}}\right)}^{\mathrm{4}}\right]u\left[{\mathit{\alpha }}_{\mathrm{m}}\left(z\right)\right]}{N_N\left(Z\right)\mathrm{\Delta }{\mathit{\sigma }}_N}}\\ & ={\displaystyle \frac{\left[\mathrm{2}-{\left(\frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{1}}}\right)}^{\mathrm{4}}-{\left(\frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{3}}}\right)}^{\mathrm{4}}\right]{\mathit{\sigma }}_{\mathrm{m}}u\left[N_{\mathrm{m}}\left(z\right)\right]}{N_N\left(Z\right)\mathrm{\Delta }{\mathit{\sigma }}_N}}\end{array}\end{array}$$

$$\begin{array}{}\text{(15)}& {\displaystyle }\begin{array}{rl}& U_{\mathrm{OAD}}\left(z\right)={\displaystyle \frac{u\left[\mathrm{OAD}\left(z\right)\right]}{N_N\left(Z\right)\mathrm{\Delta }{\mathit{\sigma }}_N}}\\ & ={\displaystyle \frac{u\left[\mathrm{2}{O}_{\mathrm{abs}}\left({\mathit{\lambda }}_{\mathrm{2}},z\right)-O_{\mathrm{abs}}\left({\mathit{\lambda }}_{\mathrm{1}},z\right)-O_{\mathrm{abs}}\left({\mathit{\lambda }}_{\mathrm{3}},z\right)\right]}{N_N\left(Z\right)\mathrm{\Delta }{\mathit{\sigma }}_N}}\\ & ={\displaystyle \frac{\left[\mathrm{2}{\mathit{\sigma }}_{\mathrm{o}}\left({\mathit{\lambda }}_{\mathrm{2}},z\right)-{\mathit{\sigma }}_{\mathrm{o}}\left({\mathit{\lambda }}_{\mathrm{1}},z\right)-{\mathit{\sigma }}_{\mathrm{o}}\left({\mathit{\lambda }}_{\mathrm{3}},z\right)\right]u\left[N_{\mathrm{o}}\left(z\right)\right]}{N_N\left(Z\right)\mathrm{\Delta }{\mathit{\sigma }}_N}}\end{array}\text{(16)}& {\displaystyle }\begin{array}{rl}{U}_{\mathrm{AED}}\left(z\right)& ={\displaystyle \frac{u\left[\mathrm{AED}\left(z\right)\right]}{N_N\left(Z\right)\mathrm{\Delta }{\mathit{\sigma }}_N}}\\ & ={\displaystyle \frac{\left[\mathrm{2}-{\left(\frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{1}}}\right)}^e-{\left(\frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{3}}}\right)}^e\right]u\left[{\mathit{\alpha }}_{\mathrm{a}}(z,s)\right]}{N_N\left(Z\right)\mathrm{\Delta }{\mathit{\sigma }}_N}}\end{array}\text{(17)}& {\displaystyle }\begin{array}{rl}& U_B\left(z\right)=\\ & {\displaystyle \frac{u\left\{\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{d}}{\mathrm{d}z}\left[\ln \frac{\begin{array}{c}\left[{\left(\frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{3}}}\right)}^{\mathrm{4}}{\mathit{\beta }}_{\mathrm{m}}\left(z\right)+{\left(\frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{3}}}\right)}^e{\mathit{\beta }}_{\mathrm{a}}(z,s)\right]\\ \times \left[{\left(\frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{1}}}\right)}^{\mathrm{4}}{\mathit{\beta }}_{\mathrm{m}}\left(Z\right)+{\left(\frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{1}}}\right)}^e{\mathit{\beta }}_{\mathrm{a}}(z,s)\right]\end{array}}{{\left[{\mathit{\beta }}_{\mathrm{m}}\left(Z\right)+{\mathit{\beta }}_{\mathrm{a}}(z,s)\right]}^{\mathrm{2}}}\right]\right\}}{N_N\left(Z\right)\mathrm{\Delta }{\mathit{\sigma }}_N}}\end{array}\end{array}$$

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& U_{\mathrm{S}}\left(z\right)={\displaystyle \frac{u\left\{\frac{\mathrm{1}}{\mathrm{2}}\times \frac{\mathrm{d}}{\mathrm{d}z}\left[\ln \frac{X\left({\mathit{\lambda }}_{\mathrm{1}},Z\right)\left({\mathit{\lambda }}_{\mathrm{3}},Z\right)}{X{\left({\mathit{\lambda }}_{\mathrm{2}},Z\right)}^{\mathrm{2}}}\right]\right\}}{N_N\left(Z\right)\mathrm{\Delta }{\mathit{\sigma }}_N}}\\ & ={\displaystyle \frac{\begin{array}{c}\frac{\mathrm{1}}{\mathrm{2}}\times \sqrt{{\left\{\frac{\mathrm{d}\left\{\frac{\mathrm{d}}{\mathrm{d}z}\left[\ln \frac{X\left({\mathit{\lambda }}_{\mathrm{1}},Z\right)\left({\mathit{\lambda }}_{\mathrm{3}},Z\right)}{X{\left({\mathit{\lambda }}_{\mathrm{2}},Z\right)}^{\mathrm{2}}}\right]\right\}}{\mathrm{d}\left[X\left({\mathit{\lambda }}_{\mathrm{1}},Z\right)\right]}\times u\left[X\left({\mathit{\lambda }}_{\mathrm{1}},Z\right)\right]\right\}}^{\mathrm{2}}}\\ \stackrel{\mathrm{‾}}{+{\left\{\frac{\mathrm{d}\left\{\frac{\mathrm{d}}{\mathrm{d}z}\left[\ln \frac{X\left({\mathit{\lambda }}_{\mathrm{1}},Z\right)\left({\mathit{\lambda }}_{\mathrm{3}},Z\right)}{X{\left({\mathit{\lambda }}_{\mathrm{2}},Z\right)}^{\mathrm{2}}}\right]\right\}}{\mathrm{d}\left[X\left({\mathit{\lambda }}_{\mathrm{2}},Z\right)\right]}\times u\left[X\left({\mathit{\lambda }}_{\mathrm{2}},Z\right)\right]\right\}}^{\mathrm{2}}}\\ \stackrel{\mathrm{‾}}{+{\left\{\frac{\mathrm{d}\left\{\frac{\mathrm{d}}{\mathrm{d}z}\left[\ln \frac{X\left({\mathit{\lambda }}_{\mathrm{1}},Z\right)\left({\mathit{\lambda }}_{\mathrm{3}},Z\right)}{X{\left({\mathit{\lambda }}_{\mathrm{2}},Z\right)}^{\mathrm{2}}}\right]\right\}}{\mathrm{d}\left[X\left({\mathit{\lambda }}_{\mathrm{3}},Z\right)\right]}\times u\left[X\left({\mathit{\lambda }}_{\mathrm{3}},Z\right)\right]\right\}}^{\mathrm{2}}}\end{array}}{N_N\left(Z\right)\mathrm{\Delta }{\mathit{\sigma }}_N}},\end{array}\end{array}$$

where $U_{{\mathrm{NO}}_{\mathrm{2}}}$ is NO₂ total retrieval relative uncertainty using the three-wavelength DIAL technique; U_MED, U_OAD, U_AED, U_B and U_S are NO₂ retrieval relative uncertainty caused by molecules, absorption of gases other than the gas of interest, aerosol (extinction and backscattering) and noise of lidar signals expressed as Eqs. (14), (15), (16), (17) and (18); u is the uncertainty function; N_m and N_o are the number density (ND) of air and ozone; σ_m is the Rayleigh scattering cross-section; σ_o is the absorption cross-section of ozone; and S is the lidar ratio.

Figure 11 NO₂ number density relative uncertainty owing to the noise of signals with an average of 1 and 2 min.

Figure 12 NO₂ number density total relative uncertainty with e = 1 (a), e = 2 (b) and e = 3 (c).

From Eqs. (14) and (15), U_MED and U_OAD are determined by u[N_m(z)] and u[N_o(z)] (uncertainties in N_m and N_o). In our measurements, profiles of temperature and pressure from local radiosondes are used to calculate N_m. Usually, one radiosonde is launched for about 8 h of measurement. One profile of air number density from local radiosondes is used to correct 8 h NO₂ measurements. According to statistics of 8 h variation in temperature and pressure in four local seasons, the uncertainty in N_m is between 1 % and 3 %. U_MED using the two-wavelength DIAL technique and the three-wavelength DIAL technique is calculated according to Eq. (14), with the uncertainty in N_a as 1 %, 2 % and 3 %, as shown in Fig. 7a. U_MED using the three-wavelength DIAL technique is far less than using the two-wavelength DIAL technique. N_o is obtained from local measurements. Because of very low values of ozone absorption cross-section differentials, with the uncertainty in N_o as 50 % and 100 %, U_OAD using the two-wavelength DIAL technique and U_OAD using the three-wavelength DIAL technique are both less than 0.5 %, as seen in Fig. 7b. Ozone absorption correction is neglected in NO₂ retrieval. From Eqs. (16) and (17), U_AED and U_B are determined by uncertainties in a_a, β_a and e. For the HU lidar system, 532 nm elastic signals are used to calculate a_a and β_a with the Fernald method to correct NO₂ retrieval (Fernald et al., 1972). A value of 50 sr is usually chosen as a lidar ratio to retrieve a_a and β_a. The lidar ratio is variable, so uncertainties in a_a and β_a are caused by the chosen lidar ratio. The range of the lidar ratio is about 30 to 70 sr for 532 nm. The uncertainty in the lidar ratio is 40 % for 50 sr. The uncertainties in a_a and β_a are calculated with uncertainty in the lidar ratio as 40 %. Finally, U_AED and U_B using the two-wavelength DIAL technique and using the three-wavelength DIAL technique are calculated with the Ångström exponent as 1, 2 and 3, as shown in Figs. 8, 9 and 10. From these figures, U_AED and U_B using the three-wavelength DIAL technique are both less than 4 %. However, U_AED below 2 km using the two-wavelength DIAL technique is more than 90 % after correction of aerosol extinction. From Eq. (18), U_S is determined by uncertainties in three-wavelength lidar signals. The uncertainties in lidar signals with an average integration time of 1 and 2 min are derived from Poisson statistics associated with the probability of detection of a repeated random event (Megie, 1985; Leblanc et al., 2016). NO₂ number density relative uncertainty owing to the noise of lidar signals with an average integration time of 1 and 2 min is obtained, as shown in Fig. 11. We can see that U_S using the two-wavelength DIAL technique is smaller than using the three-wavelength DIAL technique. With an increase in average integration time from 1 to 2 min, U_S can be effectively reduced. At last, $U_{{\mathrm{NO}}_{\mathrm{2}}}$ (the total relative uncertainties in NO₂) with e as 1, 2 and 3 is calculated, as shown in Fig. 12a, b and c.

Figure 13 HU lidar 438, 439.5 and 441 nm elastic signals measured at 21:00 LT on 13 May 2020 (a) and 22:00 LT on 27 July 2020 (c); NO₂ profiles obtained using the three-wavelength DIAL technique, two-wavelength DIAL technique and WRF-Chem model at 21:00 LT on 13 May 2020 (b) and 22:00 LT on 27 July 2020 (d).

5 Results

The three-wavelength DIAL technique was implemented by the HU lidar measurements during two cases at night, and the resulting vertical profiles are presented in Fig. 13. All NO₂ lidar measurements presented here are obtained at times with less than 10 % cloud coverage below 8 km. HU lidar 438 nm (blue line), 439.5 nm (red line) and 441 nm (black line) elastic signals measured at 21:00 LT on 13 May 2020 and 22:00 LT on 27 July 2020 are shown in Fig. 13a and c, respectively. The average integration time for these signals is 2 min. Determined from the lidar elastic signals in Fig. 13a and c, there is an existing aerosol layer between 2.2 and 3.5 km on 13 May, while 27 July presented a clean atmosphere. Figure 13b and d show retrieved NO₂ profiles using the three-wavelength DIAL technique (red line). The black error bars in Fig. 13b and d indicate the uncertainty in NO₂ retrieval calculated using Eq. (13). In Fig. 13b, the retrieved NO₂ profile between 2.2 and 3.5 km on 13 May is smooth and not affected by the aerosol layer. The NO₂ profiles (sky-blue line and purple line) were also retrieved using the conventional two-wavelength DIAL technique without and with aerosol correction, as shown in Fig. 13b, resulting in a bump between 2.2 and 3.5 km in the NO₂ profile retrieved using the two-wavelength DIAL technique. This inconsistency suggests that the two-wavelength DIAL technique cannot remove AED of the aerosol layer between 2.2 and 3.5 km, and the retrieved NO₂ profile contains AED interference. Moreover, the NO₂ retrievals below 2 km using the two-wavelength DIAL technique shown in Fig. 13b and d are more than the three-wavelength DIAL technique, suggesting that the AED of boundary aerosol was not correctly removed. Aerosol correction is very important for NO₂ retrieval using the conventional two-wavelength DIAL technique (Sasano et al., 1985). These results suggest that the proposed three-wavelength DIAL technique can effectively remove the influence of aerosol on the retrieval of NO₂. As a first-order assessment of the HU lidar NO₂profiles, we compare the retrieval results to simulated data from the Weather Research and Forecasting Chemistry (WRF-Chem) model (Grell et al., 2005) at 12 km × 12 km spatial resolution and 200 m vertical resolution. Past studies have demonstrated that WRF-Chem-simulated NO₂ results show good agreement between the OMI satellite measurements and aircraft measurements (Amnuaylojaroen et al., 2019; Barten et al., 2020), providing a data source to examine the accuracy of the HU retrievals using both the two-wavelength DIAL technique and the three-wavelength DIAL technique. The HU local NO₂ profiles for these two cases are simulated using the WRF-Chem model and shown in Fig. 13b and d. WRF-Chem-simulated NO₂ magnitudes tend to be lower compared to HU-retrieved NO₂ profiles using the three-wavelength DIAL technique (typically within ±0.1 ppb), except above 3.5 km on 13 May 2020; however, the comparison demonstrates a consistent vertical profile shape between observations and the model simulation. And retrieval results using the three-wavelength DIAL technique are much closer to simulated values compared to using the two-wavelength DIAL technique. These figures also demonstrate that the reduced fluctuations caused by aerosol backscatter when using the three-wavelength DIAL technique result in vertical profiles of NO₂ which are much more consistent with simulated data when compared to results of the two-wavelength DIAL retrievals. Both the WRF-Chem-simulated profiles and the HU retrievals of NO₂ using the three-wavelength DIAL technique are associated with uncertainties which could result in the differences in magnitude; however, the consistent nature in the vertical profile shapes from both data sources provides confidence that the HU lidar is retrieving NO₂ vertical profiles using the three-wavelength DIAL technique in the troposphere.

6 Conclusion

This study describes a lidar retrieval technique using three wavelengths simultaneously emitted from an OPO laser to measure tropospheric NO₂ profiles. The three-wavelength DIAL retrieval equations describe how the retrievals decrease errors caused by aerosol interference. Aerosol extinction differences using this proposed technique can be decreased to less than 2 % of the aerosol extinction differences resulting from a conventional two-wavelength DIAL technique. Comparing the HU lidar results to WRF-Chem model output demonstrates that the NO₂ magnitudes and vertical structure are in much better agreement with simulated data when applying the three-wavelength DIAL technique compared to using the two-wavelength technique. In the future, we will add new filters to obtain daytime NO₂ measurements. We also plan to purchase NO₂ balloon sondes for acquiring true validation data to evaluate HU lidar NO₂ results.