Retrieval of total column and surface NO₂ from Pandora zenith-sky measurements

Pandora spectrometers can retrieve nitrogen dioxide (NO₂) vertical column densities (VCDs) via two viewing geometries: direct Sun and zenith sky. The direct-Sun NO₂ VCD measurements have high quality (0.1 DU accuracy in clear-sky conditions) and do not rely on any radiative transfer model to calculate air mass factors (AMFs); however, they are not available when the Sun is obscured by clouds. To perform NO₂ measurements in cloudy conditions, a simple but robust NO₂ retrieval algorithm is developed for Pandora zenith-sky measurements. This algorithm derives empirical zenith-sky NO₂ AMFs from coincident high-quality direct-Sun NO₂ observations. Moreover, the retrieved Pandora zenith-sky NO₂ VCD data are converted to surface NO₂ concentrations with a scaling algorithm that uses chemical-transport-model predictions and satellite measurements as inputs. NO₂ VCDs and surface concentrations are retrieved from Pandora zenith-sky measurements made in Toronto, Canada, from 2015 to 2017. The retrieved Pandora zenith-sky NO₂ data (VCD and surface concentration) show good agreement with both satellite and in situ measurements. The diurnal and seasonal variations of derived Pandora zenith-sky surface NO₂ data also agree well with in situ measurements (diurnal difference within ±2 ppbv). Overall, this work shows that the new Pandora zenith-sky NO₂ products have the potential to be used in various applications such as future satellite validation in moderate cloudy scenes and air quality monitoring.

1 Introduction

Nitrogen dioxide (NO₂) is an important air pollutant and plays a critical role in tropospheric photochemistry (e.g., ECCC, 2016; EPA, 2014). It is primarily emitted from combustion processes such as fossil fuel combustion (e.g., traffic, electricity generation from power plants) and biomass burning, as well as from lightning. NO₂ is a nitrate aerosol precursor, and it also contributes to acid deposition and eutrophication (ECCC, 2016). Exposure to NO₂ can lead to adverse health effects, such as irritation of the lungs, a decrease in lung function, and an increase in susceptibility to allergens for people with asthma (EEA, 2017; WHO, 2017).

As surface NO₂ concentrations are regulated by many environmental agencies (e.g., Environment and Climate Change Canada and US Environment Protection Agency), in situ NO₂ measurements are commonly carried out by many national monitoring networks, such as the National Air Pollution Surveillance (NAPS; https://www.canada.ca/en/environment-climate-change/services/air-pollution/monitoring-networks-data/national-air-pollution-program.html, last access: 15 August 2019) network in Canada, which was established in 1969. The in situ methods used to measure surface NO₂ have evolved over the years; for example, luminol chemiluminescence (e.g., Kelly et al., 1990; Maeda et al., 1980; Wendel et al., 1983), long-path differential optical absorption spectroscopy (e.g., Platt, 1994), photolytic conversion/chemiluminescence (e.g., Gao et al., 1994; Ryerson et al., 2000), and laser-induced fluorescence (e.g., Thornton et al., 2000) are all found to be reliable methods with an uncertainty within 10 % at 1 ppbv and higher concentration levels (McClenny, 2000). Currently, the in situ approach used by NAPS for surface NO₂ air quality monitoring is the photolytic conversion/chemiluminescence technique, which converts NO₂ to NO and subsequently detects the NO by chemiluminescence reaction (McClenny, 2000; NRC, 1992). This in situ monitoring measurements provides good measurements at ground level (0.4 ppbv accuracy), but NO₂ is not uniformly mixed through the atmosphere, and not even within the atmospheric boundary layer due to emission and removal processes taking place at the surface.

Total vertical column NO₂ can be measured by many ground-based UV–visible remote-sensing instruments using direct-Sun, zenith-sky, or off-axis spectroscopy techniques (Cede et al., 2006; Drosoglou et al., 2017; Herman et al., 2009; Lee et al., 1994; Noxon, 1975; Piters et al., 2012; Roscoe et al., 2010; Tack et al., 2015; Vaughan et al., 1997). These measurements are of high quality and good precision, and have been widely used for atmospheric chemistry studies (e.g., Adams et al., 2012; Hendrick et al., 2014) and satellite validations (e.g., Celarier et al., 2008; Drosoglou et al., 2018; Irie et al., 2008; Wenig et al., 2008). Among all these different viewing geometries, direct-Sun measurements are of high accuracy and are not dependent on radiative transfer models (RTMs) to calculate air mass factors (AMFs) (Herman et al., 2009) or on knowledge of other atmospheric constituents. Zenith-sky observations have been widely used for stratospheric ozone and NO₂ observations, particularly under cloudy conditions when direct-Sun measurements are unreliable (note that zenith-sky observations use scattered sunlight and are less sensitive to clouds, e.g., Zhao et al., 2019). Off-axis measurements have good sensitivity in the boundary layer and could provide tropospheric trace gas profiles and surface concentrations (Frieß et al., 2011; Hendrick et al., 2014; Kramer et al., 2008; Wagner et al., 2011), but they are more sensitive to cloud cover than zenith-sky measurements.

The Pandora Sun spectrometer is a new instrument developed to measure vertical column densities (total columns) of trace gases in the atmosphere using Sun and sky radiation in the UV–visible part of the spectrum (Herman et al., 2009). One of its primary data products is NO₂ total vertical column density (VCD) from the direct-Sun viewing mode, where VCD represents the vertically integrated number of molecules per unit area and is reported in units of molec cm⁻² or Dobson units (1 DU $=\mathrm{2.6870}\times {\mathrm{10}}^{\mathrm{16}}$ molec cm⁻²). The Pandora direct-Sun NO₂ VCD products have been validated through many field campaigns (Flynn et al., 2014; Lamsal et al., 2017; Martins et al., 2016; Piters et al., 2012; Reed et al., 2015), ground-based comparisons (Herman et al., 2009; Wang et al., 2010), and satellite validations (Ialongo et al., 2016; Lamsal et al., 2014).

Since their introduction in 2006, Pandora spectrometers have been deployed at more than 50 sites globally. The Pandora no. 103 instrument used in this study has been deployed in Toronto, Canada since 2013 to perform direct-Sun measurements (Zhao et al., 2016). Since 2015, the observation schedule of Pandora no. 103 has been modified to perform alternating direct-Sun and zenith-sky measurements. Knepp et al. (2017) assessed Pandora's capability to derive stratospheric NO₂ using zenith-sky viewing geometry (in twilight periods), but their study was limited to slant column densities (SCDs). At this time, there are no standard Pandora zenith-sky NO₂ VCD data products available. As one goal of this work, we have focused on developing a new NO₂ retrieval algorithm for zenith-sky measurements to expand Pandora NO₂ measurements into cloudy scenes.

In addition to retrieval of zenith-sky total column NO₂, another goal of this work is to derive surface NO₂ concentration from total column measurements. Surface NO₂ has been a focus of scientific studies due to its strong correlation with air quality (AQ) and health issues (ECCC, 2016), with NO₂ as one of the three components (along with ozone and PM_2.5) used to compute the Air Quality Health Index (AQHI; Stieb et al., 2008) in Canada's AQ public awareness programs. Efforts to link total column NO₂ with its surface concentrations have been made by many researchers (Flynn et al., 2014; Knepp et al., 2015; Kollonige et al., 2017; Lamsal et al., 2008, 2014; McLinden et al., 2014). For example, Knepp et al. (2015) proposed a method to estimate NO₂ surface mixing ratios from Pandora direct-Sun total column NO₂ via application of a planetary boundary layer (PBL) height correction factor. Kollonige et al. (2017) adapted this method and compared Pandora direct-Sun surface NO₂ and Ozone Monitoring Instrument (OMI) surface NO₂. They concluded that the two main sources of error for the conversion of the total column NO₂ to surface NO₂ are (1) poor weather conditions (e.g., cloud cover and precipitation) and (2) PBL height estimation, both of which affect the NO₂ column–surface relationship and instrument sensitivities to boundary layer NO₂. Thus, in this work, we present a simple but robust algorithm for deriving surface NO₂ concentration from Pandora zenith-sky measurements, which has several advantages, such as the ability (1) to extend Pandora NO₂ measurements to cloudy conditions and (2) to provide more accurate surface NO₂ concentration estimates that are less sensitive to PBL height. This work also provides reliable total column NO₂ measurements in cloudy conditions and could be used in satellite validations in partially cloudy scenes.

This paper is organized as follows. Section 2 describes the measured and modelled NO₂ data used in this study. In Sect. 3, the empirical AMFs for Pandora zenith-sky NO₂ measurements are derived using high-quality Pandora direct-Sun total column NO₂ data. These empirical AMFs and the Network for the Detection of Atmospheric Composition Change (NDACC) AMFs (Hendrick et al., 2011; Sarkissian et al., 1995; Van Roozendael et al., 1998; Van Roozendael and Hendrick, 2009; Vaughan et al., 1997) are both applied to Pandora zenith-sky total column NO₂ retrievals to help evaluate the performance of the empirical AMFs. Also, the retrieved Pandora zenith-sky total column NO₂ data are evaluated by comparison with satellite measurements. In Sect. 4, the zenith-sky total column NO₂ data are converted to surface concentration by using a scaling algorithm. The zenith-sky surface NO₂ concentration data are assessed by comparison with in situ measurements. Lastly, in Sect. 5, several aspects of this zenith-sky surface NO₂ dataset are discussed, which include diurnal and seasonal variation, and PBL effect, followed by conclusions in Sect. 6.

2 Datasets and models

2.1 Measurements

2.1.1 Pandora direct-Sun total column NO₂

The Pandora instrument records spectra between 280 and 530 nm with a resolution of 0.6 nm (Herman et al., 2009, 2015; Tzortziou et al., 2012). It uses a temperature-stabilized Czerny–Turner spectrometer, with a 50 µm entrance slit, 1200 groove mm⁻¹ grating, and a 2048×64 back-thinned Hamamatsu charge-coupled device (CCD) detector. The spectra are analyzed using a total optical absorption spectroscopy (TOAS) technique (Cede, 2019), in which absorption cross sections for multiple atmospheric absorbers, such as ozone, NO₂, and sulfur dioxide (SO₂), are fitted to the spectra.

The Pandora direct-Sun total column NO₂ data are produced using Pandora's standard NO₂ algorithm implemented in the BlickP software (Cede, 2019). The measured direct-Sun spectra from 400 to 440 nm are used in the TOAS analysis. A synthetic reference spectrum is produced by averaging multiple measured spectra and corrected for the estimated total optical depth included in it. Cross sections of NO₂ at an effective temperature of 254.5 K (Vandaele et al., 1998), ozone at an effective temperature of 225 K (Brion et al., 1993, 1998; Daumont et al., 1992), and a fourth-order polynomial are all fitted. The resulting NO₂ SCDs are then converted to total column VCDs by using direct-Sun geometry AMFs. Herman et al. (2009) show that Pandora direct-Sun total column NO₂ has a clear-sky precision of 0.01 DU (in the slant column) and a nominal accuracy of 0.1 DU (in the vertical column, 2σ level). Additional information on Pandora calibrations, operation, and retrieval algorithms can be found in Herman et al. (2009) and Cede (2019).

The Pandora no. 103 instrument has been deployed in Toronto since September 2013 to perform direct-Sun observations (Zhao et al., 2016). The instrument is installed on the roof of the Environment and Climate Change Canada (ECCC) Downsview building (43.7810^∘ N, −79.4680^∘ W) in Toronto. The building is located in a suburban area with multiple roads nearby. Since 2015, the instrument has been employing an alternating direct-Sun and zenith-sky observation schedule, which consists of direct-Sun measurements every 90 s and zenith-sky measurements every 30 min during the sunlit period. About 2.5 years (February 2015 to September 2017) of continuous alternating measurements are used in this study.

2.1.2 Pandora zenith-sky total column NO₂

Retrieval of trace gases from Pandora's zenith-sky measurements is not included in the standard BlickP processing software (Cede, 2019). The Pandora zenith-sky spectra for this study are processed using the differential optical absorption spectroscopy (DOAS) technique (Noxon, 1975; Platt, 1994; Platt and Stutz, 2008; Solomon et al., 1987) with the QDOAS software (Danckaert et al., 2015). A single reference spectrum is used, which was obtained from a zenith-sky measurement at local noon from a day that had low total column NO₂. Following the NDACC recommendations (Van Roozendael and Hendrick, 2012), NO₂ differential slant column densities (dSCDs) are retrieved in the 425–490 nm window (to retrieve oxygen collision complex simultaneously). The oxygen collision complex (O₂)₂ (referred here as O₄), which is created by the collision of two oxygen molecules, has broadband absorptions from UV to near-IR spectral ranges (Greenblatt et al., 1990; Platt and Stutz, 2008; Thalman and Volkamer, 2013). O₄ is widely used as a reference gas by many DOAS applications to infer cloud and aerosol properties (e.g., Gielen et al., 2014; Wagner et al., 2004, 2014, 2016, 2019; Wang et al., 2015; Zhao et al., 2019). Cross sections of NO₂ at an effective temperature of 254.5 K (Vandaele et al., 1998), ozone at an effective temperature of 223 K (Bogumil et al., 2003), H₂O (Rothman et al., 2005), O₄ (Hermans et al., 2003), and ring (Chance and Spurr, 1997) are all fitted; a fifth-order polynomial and a first-order linear offset are also included in the DOAS analysis.

The output of QDOAS is NO₂ dSCDs, which can be converted to total column NO₂ via the Langley plot method with the use of the NDACC NO₂ AMF look-up table (LUT) (Van Roozendael and Hendrick, 2012). The NDACC AMF LUT is used here only as a reference since it was primarily developed for retrieval of stratospheric NO₂. Other empirical zenith-sky NO₂ AMFs have been developed and are used to convert NO₂ dSCDs to total columns. Details about these two different AMFs are given in Sect. 3.1.

2.1.3 OMI SPv3 data

OMI is a Dutch–Finnish nadir-viewing UV–visible spectrometer aboard the National Aeronautics and Space Administration (NASA)'s Earth Observing System (EOS) Aura satellite that was launched in July 2004. The OMI instrument measures the solar radiation backscattered by the Earth's atmosphere and surface between 270 and 500 nm with resolution of 0.5 nm (Levelt et al., 2006, 2018). OMI has a 780×576 CCD detector that measures at 60 across-track positions simultaneously and thus does not require across-track scanning. Due to this approach, the spatial resolution of the CCD pixels varies significantly along the across-track direction: those pixels near the track centre have a ground footprint of 13 km×24 km (along track × across track), whereas those close to the track edge (e.g., view zenith angle =56^∘) have a ground footprint roughly of 23 km×126 km (de Graaf et al., 2016). Note that from 2012 onwards the smallest pixels (across-track positions) can no longer be used and are excluded from the analysis (known as the “row anomaly”, i.e., Levelt et al., 2018). This means the “smallest” pixels available for an OMI comparison are larger than 13 km×24 km.

The OMI NO₂ data used in this work are the NASA standard product (SP) (Bucsela et al., 2013; Wenig et al., 2008) version 3.0 level 2 (SPv3.0) (Krotkov et al., 2017). The NO₂ SCDs are derived using the DOAS technique in the 405–465 nm window (Marchenko et al., 2015). The AMFs used in SPv3.0 are calculated by using $\mathrm{1}{}^{\circ }\times \mathrm{1.25}{}^{\circ }$ (latitude × longitude) resolution a priori NO₂ and temperature profiles from the Global Modeling Initiative (GMI) chemistry–transport model with yearly varying emissions (Krotkov et al., 2017).

2.1.4 In situ measurements

The NAPS network was established in 1969 to monitor and assess the quality of ambient (outdoor) air in the populated regions of Canada. NAPS provides accurate long-term air quality data (ozone, NO₂, SO₂, carbon monoxide (CO), fine particulate matter, etc.) of a uniform standard across Canada (e.g., Dabek-Zlotorzynska et al., 2011; Reid and Aherne, 2016).

The in situ NO₂ data used in this study were collected at the NAPS Toronto north station (located 100 m away from the Pandora instrument). The site is 186 m above sea level, and the height of the air intake is 4 m above the ground.

The in situ NO₂ concentration is measured using a photolytic NO₂ instrument (Thermo 42i) that is also sensitive to other gaseous inorganic nitrogen compounds (e.g., nitric acid (HNO₃) and peroxyacetyl nitrate (PAN)) (McLinden et al., 2014). Thus, in areas where direct NO_x (nitrogen oxides) emission sources are limited and other nitrogen compounds are present, NO₂ may be overestimated (e.g., in rural areas). For the current site, however, this positive bias has been found to be only about 5 %, except for very low NO₂ concentrations (< 5 ppbv) (Yushan Su, Ontario Ministry of the Environment, Conservation and Parks, personal communication, October 2018).

2.2 Numerical models

Predicted NO₂ fields from three atmospheric chemistry models are used in the algorithm described in Sect. 4.1 to derive surface NO₂ concentration from Pandora zenith-sky total column NO₂ data. Following McLinden et al. (2014), this work uses the Global Environmental Multi-scale Modelling Air quality and CHemistry (GEM-MACH) regional chemical transport model (CTM) and the GEOS-Chem global CTM to simulate total columns and vertical profiles of tropospheric NO₂ and surface NO₂ concentration. The stratospheric NO₂ partial columns are estimated using OMI satellite data and the Pratmo box model.

2.2.1 GEM-MACH

GEM-MACH is ECCC's regional air quality forecast model. It is run operationally twice per day to predict hourly surface pollutant concentrations over North America for the next 48 h (Moran et al., 2009; Pavlovic et al., 2016; Pendlebury et al., 2018). The model consists of an online tropospheric chemistry module (Akingunola et al., 2018; Pavlovic et al., 2016) embedded within the ECCC Global Environmental Multi-scale (GEM) numerical weather prediction model (Côté et al., 1998). Physical and chemical processes represented in GEM-MACH include emissions, dispersion, gas- and aqueous-phase chemistry, inorganic heterogeneous chemistry, aerosol dynamics, and wet and dry removal. The model uses gridded hourly emission fields based on US and Mexican national inventories from the US Environmental Protection Agency (EPA) Air Emissions Modeling Platform and on Canada's national Air Pollutant Emission Inventory (APEI; https://pollution-waste.canada.ca/air-emission-inventory, last access: 23 November 2018) (Zhang et al., 2018). Currently, only NO_x emissions in the PBL are included in the operational model; free-tropospheric NO_x emissions from lightning and in-flight aircraft are not considered. In this work, the GEM-MACH hourly NO₂ vertical profiles from 0 to 1.5 km and surface concentrations are retrieved from archived operational forecasts on the native model grid covering North America at 10 km×10 km horizontal resolution for the period April 2016 to December 2017. The corresponding grid box closest to the Pandora location was used in this study.

2.2.2 GEOS-Chem

The GEOS-Chem chemical transport model (Bey et al., 2001) has been used extensively in the retrieval of tropospheric columns and has been shown to be capable of reasonably simulating the vertical distributions of NO₂ (Lamsal et al., 2008; Martin et al., 2002; McLinden et al., 2014). The model has a detailed representation of tropospheric chemistry, including aerosols and their precursors (Park et al., 2004). In the simulation used in this study, a global lightning NO_x source of 6 Tg N yr⁻¹ (Martin et al., 2002) was imposed. Lightning NO_x emissions are computed as a function of cloud-top height and are scaled globally as described by Sauvage et al. (2007) to match Optical Transient Detector/Lightning Imaging Sensor (OTD/LIS) climatological observations of lightning flashes. The model was run on a $\mathrm{1}/\mathrm{2}{}^{\circ }\times \mathrm{2}/\mathrm{3}{}^{\circ }$ (latitude × longitude) grid in nested mode over North America and was driven by assimilated meteorology from the Goddard Earth Observing System (GEOS-5). The modelled NO₂ profiles were used to calculate monthly mean NO₂ partial columns in the free troposphere (1.5 to 12 km), as the GEM-MACH model does not include free-tropospheric NO₂ sources (lightning, in-flight aircraft emissions).

2.2.3 Pratmo box model

Pratmo is a stratospheric photochemical box model (Brohede et al., 2008; Lindenmaier et al., 2011; McLinden et al., 2000). The model has detailed stratospheric chemistry that includes long-lived species (nitrous oxide (N₂O), methane (CH₄), and water vapor (H₂O)) and halogen families (NO_y, Cl_y, and Br_y) that are based on a combination of three-dimensional model output and tracer correlations (Adams et al., 2017). Heterogeneous chemistry of background stratospheric sulfate aerosols is also included. The model is constrained with climatological profiles of ozone and temperature.

Stratospheric NO₂ has a strong diurnal variation; therefore, diurnal corrections must be applied when OMI stratospheric NO₂ measurements (around local noon) are interpolated to Pandora measurement times. Ratios of modelled stratospheric NO₂ columns are calculated at OMI overpass time and Pandora measurement time. These ratios are multiplied by the OMI measured stratospheric NO₂ to produce stratospheric NO₂ columns corresponding to the time of Pandora measurements. Details about the use of the Pratmo box model and the calculation of stratospheric NO₂ partial columns are provided in Sect. 4.1.

3 Total column NO₂ retrieval

3.1 Zenith-sky air mass factor

The NDACC UV–visible network uses zenith-sky AMFs in its total column NO₂ retrievals. To improve the overall homogeneity of the UV–visible NO₂ column measurements, NDACC recommended using the NO₂ AMF LUT (Van Roozendael and Hendrick, 2012). This LUT is based on climatological NO₂ profiles that are composed of (1) 20–60 km NO₂ profiles developed by Lambert et al. (1999, 2000) and (2) 12–20 km NO₂ profiles derived from SAOZ (Système D'Analyse par Observations Zénithales) balloon observations (Van Roozendael and Hendrick, 2012). The NO₂ concentration is set to zero below 12 km altitude. The NO₂ AMFs have been calculated using the UVSPEC/DISORT RTM (Hendrick et al., 2006; Wagner et al., 2007). The parameters used in building the LUT are wavelength, ground albedo, altitude of the station, and solar zenith angle (SZA). Aerosol extinction, ozone, and temperature profiles come from an aerosol model (Shettle, 1989), the US Standard Atmosphere, and the TOMS V8 climatology, respectively.

The NDACC LUT is designed for stratospheric NO₂ retrievals. Note that the absence of tropospheric NO₂ in the NDACC LUT construction will lead to an underestimation of the total column NO₂ in urban areas. For example, from 2015 to 2017, tropospheric NO₂ accounted for 73±11 % (1σ) of the total column amounts in Toronto (OMI SPv3.0 data). To account for this significant tropospheric NO₂ in urban areas, new empirical AMFs were developed in this study and the NDACC AMF LUT is used for comparison purposes only. In Tack et al. (2015), a more sophisticated four-step approach to derive total and tropospheric NO₂ columns from zenith-sky measurements was proposed, which involved using a RTM to calculate appropriate tropospheric AMFs. However, due to benefits from using the high-quality Pandora direct-Sun total column NO₂ measurements, this work took a different but simple and robust approach to derive zenith-sky total column NO₂.

Empirical AMFs are calculated for Pandora zenith-sky NO₂ measurements in such a way that they can be used to retrieve zenith-sky total column NO₂ values that match the high-quality Pandora direct-Sun total column NO₂ values. Inferring total columns from zenith-sky observations through comparisons with accurate direct-Sun observations is a common approach for Brewer and Dobson zenith-sky total ozone measurements (Kerr et al., 1988). For example, in the Brewer instrument zenith-sky ozone algorithm, weighted zenith-sky light intensities measured at four wavelengths (F) are expressed as a function of the slant path (μ) and total column ozone (Kerr et al., 1981). The nine semi-empirical coefficients used to derive total column ozone from measured F in the equation are estimated from a set of direct-Sun and zenith-sky observations made nearly simultaneously (Fioletov et al., 2011). Instead of finding the link between zenith-sky spectral intensity and total column values (i.e., following the Brewer and Dobson zenith-sky total ozone retrieval method), deriving empirical zenith-sky AMFs for Pandora zenith-sky measurements is more straightforward since Pandora zenith-sky spectra can be analyzed to produce NO₂ dSCDs.

The relation between VCD and dSCD can be expressed as

$$\begin{array}{}\text{(1)}& \mathrm{VCD}={\displaystyle \frac{\mathrm{dSCD}+\mathrm{RCD}}{\mathrm{AMF}}},\end{array}$$

where RCD is the reference column density that shows the slant column amount of the trace gas in the reference spectrum (Sect. 2.1.2). If we make an assumption that the coincident direct-Sun (DS) and zenith-sky (ZS) measurements sampled the same air mass, then the empirical zenith-sky AMFs (referred to here as AMF_ZS-Emp) can be calculated by assuming VCD_DS= VCD_ZS, which gives

$$\begin{array}{}\text{(2)}& {\mathrm{VCD}}_{\mathrm{DS}}\left(\mathrm{SZA}\right)={\displaystyle \frac{{\mathrm{dSCD}}_{\mathrm{ZS}}\left(\mathrm{SZA}\right)+{\mathrm{RCD}}_{\mathrm{ZS}}}{{\mathrm{AMF}}_{\mathrm{ZS}-\mathrm{Emp}}\left(\mathrm{SZA}\right)}}.\end{array}$$

Next, we can use nearly coincident VCD_DS and dSCD_ZS in a multi-non-linear regression to retrieve AMF_ZS-Emp and RCD_ZS together. To ensure the quality of the retrieved AMF_ZS-Emp, only high-quality direct-Sun total column NO₂ data are used with SZA < 75^∘. Details about the empirical zenith-sky AMF calculation are shown in Appendix A.

Figure 1 shows a comparison of the empirical zenith-sky AMFs and NDACC AMFs (calculated for the Toronto measurements). Total column NO₂ can then be retrieved using Eq. (1) and these two sets of AMFs, where the one based on empirical AMFs is referred to as VCD_ZS-Emp and the one based on NDACC AMFs is referred to as VCD_ZS-NDACC. The RCD value used in the retrievals is 0.39±0.01 DU, which is retrieved along with AMF_ZS-Emp (Appendix A). Figure 2 shows the comparisons of the NO₂ columns measured by zenith-sky and direct-Sun methods. The regression analyses were performed by using the following coincidence criteria: (1) nearest Pandora direct-Sun measurement that was within ±5 min of Pandora zenith-sky measurement, (2) SZA < 75^∘, and (3) Pandora direct-Sun total column NO₂ data have assured high quality (BlickP L2 data quality flag for nitrogen dioxide is 0). In general, the VCD_ZS-Emp and VCD_ZS-NDACC performed as expected. Compared with VCD_DS, the VCD_ZS-NDACC shows a −25 % bias, while the VCD_ZS-Emp only shows a −4 % bias (indicated by the red lines on each panel and their slopes). In addition, VCD_ZS-Emp shows less SZA dependence than VCD_ZS-NDACC (see the increased bias for measurements made in larger SZA conditions in Fig. 2b). These results confirm that, for urban sites, the tropospheric NO₂ profile should be included when calculating empirical zenith-sky AMFs. In the rest of the paper, only the zenith-sky NO₂ retrieved using empirical AMFs will be discussed. The derived zenith-sky total column NO₂ values are affected by both clouds and aerosols due to their impact on the light path. The presence of clouds and aerosols contributes to the uncertainty of the measurements. However, the impact of aerosols is expected to be moderate in most cases compared to that of clouds (e.g., Hendrick et al., 2011; Tack et al., 2015). Thus, this work has focused on evaluating the impact from clouds. Note that the Pandora zenith-sky total column NO₂ data discussed in Sect. 3 are a “clear-sky subset” of Pandora zenith-sky measurements. The assessment of Pandora zenith-sky NO₂ measurements in cloudy conditions is provided in Sect. 4.

Figure 1 Comparison of zenith-sky NO₂ air mass factors. Blue and red squares with error bars (standard error) represent the empirical discrete zenith-sky NO₂ AMFs in each SZA bin for Toronto for the period February 2015 to September 2017. Blue and red lines show the fitted empirical zenith-sky NO₂ AMFs. NDACC AMFs calculated using the NDACC look-up table and assuming no NO₂ in the troposphere are shown in yellow.

Figure 2 Comparisons of NO₂ total columns (2015–2017): (a) zenith-sky total column NO₂ retrieved using empirical AMFs vs. direct-Sun total column NO₂, (b) zenith-sky total column NO₂ retrieved using NDACC AMFs vs. direct-Sun total column NO₂. On each scatter plot, the red line is the linear fit with intercept set to 0, and the black line is the one-to-one line. The scatter plot is colour-coded by solar zenith angle (SZA).

3.2 Comparison with satellite measurements

To illustrate the NO₂ variability over Toronto, Fig. 3 shows the time series (2015–2017) from Pandora direct-Sun, zenith-sky, and OMI SPv3.0 total column NO₂. In general, the NO₂ datasets from the ground-based Pandora instrument and the satellite follow the same pattern. However, the satellite data are likely to miss the peak NO₂ values in the morning since OMI only passes over Toronto once per day around 13:30 LT (local time).

Figure 3 Annual time series of Pandora direct-Sun (DS), Pandora zenith-sky (ZS), and OMI SPv3 total column NO₂ in Toronto from 2015 to 2017.

We also performed regression analyses by using the following coincidence criteria: (1) nearest (in time) measurement that was within ±30 min of OMI overpass time, (2) closest OMI ground pixel (having a distance from the ground pixel centre to the location of the Pandora instrument less than 20 km), and (3) cloud fraction < =0.3 (the effective geometric cloud fraction, as determined by the OMCLDO2 algorithm; Celarier et al., 2016). In this comparison, only high-quality OMI data are used (VcdQualityFlags =0) (Celarier et al., 2016). Figure 4a and b show the scatter plots of OMI vs. Pandora direct-Sun and OMI vs. Pandora zenith-sky total column NO₂, respectively. Figure 4c and d show similar comparisons but only use OMI NO₂ measured by “small pixels” (i.e., having viewing zenith angle of less than 35^∘). The better correlation and lower bias for zenith-sky vs. direct-Sun measurements might be a case of coincident errors; i.e., compared to Pandora direct-Sun total column NO₂, both OMI and Pandora zenith-sky total column NO₂ underestimate the local NO₂ at Toronto (see Fig. 2). When taking into account the standard error of the fitting and the confidence level of R, the difference between zenith-sky and direct-Sun data is not significant (i.e., in Fig. 4 from panels a to d, the slopes with standard error are 0.64±0.02, 0.67±0.02, 0.70±0.04, and 0.71±0.03; the 95 % confidence intervals for R values are 0.45 to 0.63, 0.61 to 0.75, 0.43 to 0.77, and 0.60 to 0.86). The comparison results indicate that, at the Toronto site, OMI underestimates the total column by about 30 %. This underestimation is qualitatively consistent with the fact that the Pandora location is near the northern edge of peak Toronto NO₂, and the relatively large OMI pixels are also generally sampling areas of less NO₂ in the vicinity. The use of the relatively coarse (1^∘) GMI model for profile shapes (Sect. 2.1.3) will also lead to a low bias considering the peak NO_x emissions span roughly $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$. Similar results have been found elsewhere.

Figure 4 OMI vs. Pandora total column NO₂ (2015–2017). Panels (a) and (c) show OMI vs. Pandora direct-Sun NO₂, and (b) and (d) show OMI vs. Pandora zenith-sky NO₂. Panels (a) and (b) show all available OMI measurements, while panels (c) and (d) show OMI data from small pixels only. On each scatter plot, the red line is the linear fit with intercept set to 0 and the black line is the one-to-one line. All scatter plots are colour-coded by the distance from the centre of an OMI ground pixel to the location of Pandora.

Ialongo et al. (2016) reported a similar negative bias using OMI SPv3.0 and Pandora direct-Sun total column NO₂ in Helsinki (−32 % bias and R=0.51), and they suggested this was due to the difference between the OMI pixel and the relatively small Pandora field of view. In Reed et al. (2015), Pandora measurements at 11 sites were evaluated; the authors found that the best correlation between OMI SPv3.0 and Pandora direct-Sun total column NO₂ data is for rural sites. They concluded this could be due to smaller atmospheric variability in the rural region. Other studies such as Goldberg et al. (2017) found an even worse OMI–Pandora comparison between these two data products with striking negative bias at high values and poor correlation (R=0.3). The authors attributed the poor agreement to the coarse resolution of OMI and its AMFs computed with GMI a priori NO₂ profiles. In general, our comparison results show that (1) the Pandora direct-Sun total column NO₂ data measured in Toronto have a reasonable agreement with OMI, and (2) the Pandora zenith-sky total column NO₂ data show results similar to those for direct-Sun total column when compared with OMI SPv3.0.

4 Surface NO₂ concentration retrieval

The performance of the clear-sky Pandora zenith-sky total column NO₂ data has been assessed by using OMI and Pandora direct-Sun data as described in Sect. 3.2. However, the validation of cloudy-scene Pandora zenith-sky total column data is not simple, since near-simultaneous good-quality direct-Sun or satellite measurements in most cloudy conditions are not available. This cloudy-scene validation can be done by comparison with in situ NO₂ measurements that are not affected by weather. In general, the comparison between total column and surface concentrations can be done by two approaches: (1) convert Pandora zenith-sky total columns to surface concentrations; and (2) convert in situ surface concentrations to total column values. For example, Spinei et al. (2018) calculated “ground-up” VCDs from in situ surface concentrations by using additional measurements of PBL height or assuming trace gas profiles. In this work, the first approach is employed since the surface NO₂ data products from Pandora remote-sensing measurements have direct applications in areas such as air quality monitoring.

4.1 Column-to-surface conversion algorithm

A simple but robust scaling method is adapted to derive surface NO₂ concentration from Pandora zenith-sky total column NO₂ measurements. Following Lamsal et al. (2008) and McLinden et al. (2014), the surface NO₂ concentration is estimated using the modelled profile and surface concentration:

$$\begin{array}{}\text{(3)}& C_{\mathrm{pan}}=\left(V_{\mathrm{pan}}-V_{\mathrm{strat}}-V_{\mathrm{ftrop}}\right)\times {\left({\displaystyle \frac{C}{V_{\mathrm{PBL}}}}\right)}_{\text{G-M}},\end{array}$$

where C_pan is the surface NO₂ volume mixing ratio (VMR) to be estimated, C is the surface NO₂ VMR from GEM-MACH (or G-M), V_pan is the total column NO₂ measured by Pandora, V_strat is the stratospheric NO₂ partial column, V_ftrop is the NO₂ partial column in the free troposphere, and V_PBL is the NO₂ partial column in the PBL. This equation assumes the chemical transport models can effectively capture the spatial and temporal behaviour of the concentration-to-partial-column ratio.

In this work, V_PBL (0–1.5 km) is integrated from the GEM-MACH NO₂ profile and V_ftrop (1.5–12 km) is integrated from the GEOS-Chem NO₂ profile. Both GEM-MACH and GEOS-Chem have an hourly temporal resolution. Thus, the integrated V_PBL and V_ftrop can account for NO₂ diurnal variation. However, V_strat is from OMI monthly mean stratospheric NO₂, which does not have diurnal variation. Thus, the Pratmo box model is used to calculate stratospheric NO₂ diurnal ratios. The OMI stratospheric NO₂ columns are interpolated to morning and evening hours by multiplying by the box-model diurnal ratios. Details about the calculation of V_strat as well as references are provided in Appendix B.

The (C∕V_PBL)_G-M ratio in Eq. (3) is provided by GEM-MACH, and has hourly temporal resolution. This modelled (C∕V_PBL)_G-M ratio is referred to here as a conversion ratio R_CV. Besides the hourly modelled conversion ratio, a simple monthly look-up table is built using an average of the 1.5 years of GEM-MACH model outputs (April 2016 to December 2017) that were available. The look-up table (referred to here as the Pandora surface-concentration look-up table, or PSC-LUT) is composed of monthly conversion ratios with hourly resolution as shown in Fig. 5. For example, assuming that a Pandora NO₂ total column measurement is made on a day in December at 15:00 LST, then the corresponding conversion ratio from the PSC-LUT is 28 ppbv DU⁻¹ (see the black arrow). Our results in Fig. 5 show that the conversion ratio changes throughout the day as well as with season: 0.1 DU (partial column NO₂ in the PBL) corresponds to 5–8 pptv of surface NO₂ in the morning (08:00 LST), 2–3 pptv around local noon (13:00 LST), and 2–4 pptv in the evening (18:00 LST). In general, the variation of conversion ratios demonstrates that the surface NO₂ concentration is controlled not only by PBL height but also by both boundary layer dynamics and photochemistry. The surface NO₂ derived using the hourly modelled R_CV ratio is referred to here as C_pan-model, while the surface NO₂ derived using the monthly mean PSC-LUT is referred to here as C_pan-LUT. In general, C_pan-model is a data product that depends on daily model outputs, but C_pan-LUT only needs the pre-calculated PSC-LUT and is thus less dependent on the model. In general, the look-up table approach (C_pan-LUT) is aiming for a quick and near-real-time data delivery. Thus, to minimize year-to-year variation (e.g., from changing meteorological conditions or changing local emission patterns), for a given year, we recommend using a mean PSC-LUT that is calculated from model simulations of previous years. On the other hand, the C_pan-model is the offline, high-quality, year-specific data product that will be delivered for air quality research and other applications. Details of these two different surface NO₂ data products are discussed in the next section.

Figure 5 Dependence of the Pandora surface NO₂ concentration look-up table (PSC-LUT) on month of year and hour of day. The PSC-LUT is constructed using the GEM-MACH modelled NO₂ conversion ratios. Solid lines are monthly mean conversion ratios colour-coded by month. The shaded envelopes are the standard error of the mean.

Figure 6 Modelled and Pandora zenith-sky surface NO₂ vs. in situ NO₂ (2016–2017). Panel (a) shows the GEM-MACH modelled surface NO₂ data vs. in situ NO₂; panels (b) and (c) show the Pandora ZS surface NO₂ data vs. in situ NO₂. The Pandora ZS surface NO₂ data in panels (b) and (c) are derived using the hourly modelled conversion ratio and the monthly PSC-LUT, respectively. Panels (d) to (f) are histograms corresponding to the data in panels (a) to (c). On each scatter plot, the red line is the linear fit with intercept set to 0 and the black line is the one-to-one line. The scatter plots are colour-coded by the normalized density of the points.

4.2 Comparison with measurements and model

Figure 6 shows the evaluation of modelled and Pandora zenith-sky surface NO₂ concentrations, both using in situ NO₂ measurements as the reference. The Pandora data have been filtered for heavy clouds (details are given in Sect. 4.3). The GEM-MACH modelled surface concentrations in Toronto reproduce the in situ measurements very well with the comparison showing high correlation (R=0.78) and moderate positive bias (37 %, Fig. 6a). The Pandora zenith-sky surface NO₂ data, C_pan-model, shows almost the same correlation (R=0.77), with only −7 % bias (Fig. 6b). The better performance of C_pan-model is expected since the conversion method for Pandora zenith-sky measurements relies on the GEM-MACH modelled NO₂ profile (see Eq. 3); in other words, the Pandora zenith-sky surface NO₂ has at least one more piece of information (i.e., NO₂ total column) than GEM-MACH surface NO₂ concentrations. The C_pan-LUT shows a similar correlation coefficient (R=0.73) and has improved bias (−3 %, Fig. 6c). This result (slightly lower correlation) is also reasonable and acceptable since C_pan-LUT is derived with the monthly PSC-LUT, which has less accurate information than the hourly modelled data.

Besides the improved bias, Pandora zenith-sky surface NO₂ concentrations, C_pan-model and C_pan-LUT (Fig. 6e and f), also have better frequency distributions than the GEM-MACH (Fig. 6d). Figure 6d shows that the NO₂ surface concentrations peaks (ambient background concentrations) from model and in situ data are misaligned. This indicates that the GEM-MACH NO₂ background surface concentrations have a 1 ppbv low bias at this site. In contrast, the zenith-sky surface NO₂ at peak frequency matches the in situ data (Fig. 6e and f), indicating that the low bias of the background surface NO₂ value has been corrected with this additional information from Pandora zenith-sky total column measurements. In addition, in high NO₂ concentration conditions (> 20 ppbv), the zenith-sky surface NO₂ also shows better agreement with the in situ NO₂ than do the modelled data. The mean of the top 10 % of the in situ data is 26±1 ppbv (uncertainty of the mean), whereas the corresponding values for GEM-MACH, C_pan-model, and C_pan-LUT are 39±1, 26±1, and 27±1 ppbv, respectively.

Figure 7 Modelled and Pandora direct-Sun surface NO₂ vs. in situ NO₂ (2016–2017). Panel (a) shows the GEM-MACH modelled surface NO₂ data vs. in situ NO₂; panels (b) and (c) show the Pandora DS surface NO₂ data vs. in situ NO₂. The Pandora DS surface NO₂ data in panels (b) and (c) are derived using the hourly modelled conversion ratio and the monthly PSC-LUT, respectively. Panels (d) to (f) are histograms corresponding to the data in panels (a) to (c). On each scatter plot, the red line is the linear fit with intercept set to 0 and the black line is the one-to-one line. The scatter plots are colour-coded by the normalized density of the points.

The total column-to-surface concentration conversion algorithm has also been applied to the Pandora direct-Sun total column NO₂ (see Fig. 7). Figure 7b shows that the direct-Sun surface NO₂ data have a similar agreement with the in situ data (−8 % bias and R=0.80) as the zenith-sky surface NO₂. In high NO₂ concentration conditions, direct-Sun data have a similarly good agreement with the in situ measurements. For this direct-Sun based dataset, the mean of the top 10 % of the in situ data is 27±1 ppbv, whereas the corresponding values for GEM-MACH, C_pan-model, and C_pan-LUT are 40±1, 27±1, and 27±1 ppbv, respectively.

Thus, in general, both Pandora zenith-sky and direct-Sun surface NO₂ datasets can be used reliably to obtain surface concentrations. The good consistency between C_pan-model and C_pan-LUT implies that two versions of Pandora surface NO₂ data can be delivered in the future, i.e., an offline version that relies on the inputs from hourly model and a near-real-time version that only needs a pre-calculated LUT.

4.3 Measurements in different sky conditions

Although zenith-sky observations are less sensitive to cloud conditions than direct-Sun observations, we still need to be cautious about the derived zenith-sky surface NO₂ in heavy cloud conditions. Due to enhanced scattering, heavy clouds could lead to a significant overestimation of surface NO₂ derived from zenith-sky measurements. A cloud filtering method based on retrieved O₄ dSCDs is used to identify these conditions. High retrieved O₄ values correspond to long optical path lengths, and therefore it is expected that corresponding NO₂ values are overestimated as discussed in Appendix C.

Figure 8 Example of surface NO₂ concentration time series in all conditions (April 2017). The in situ, Pandora DS, and Pandora ZS surface NO₂ concentrations are shown by different coloured dots. The TSI relative strength of direct-Sun data is plotted as a colour-coded horizontal dotted line in the top area of each panel. For Pandora zenith-sky data, the measurements with enhanced O₄ (heavy cloud indicator) are also labelled by green squares. Dates are in mm/dd format.

The effectiveness of the zenith-sky NO₂ in cloudy scenes is demonstrated by the time series plots (Fig. 8) of in situ and Pandora direct-Sun and zenith-sky data (in their original temporal resolutions). Under clear-sky conditions (for example, 8–14 April), both Pandora direct-Sun and zenith-sky-based surface concentrations correlate well with the in situ measurements. Under moderately cloudy conditions, when Pandora direct-Sun observations cannot provide high-quality data, Pandora zenith-sky observation still can yield good measurements that compare well with in situ data (for example, 26–29 April). Under heavy cloud conditions, however, which are identified by enhanced O₄ (Appendix C), Pandora zenith-sky-derived surface NO₂ yielded higher than in situ measurements (for example, 4 and 6 April; see the green squares). This feature is due to the enhanced multi-scattering in heavy cloud conditions, which leads to enhanced NO₂ absorption in the measured spectra.

Sensitivity tests (Appendix C) show that only 10 % of all zenith-sky measurements are strongly affected by this enhanced absorption, indicating the zenith-sky NO₂ algorithm is applicable to most measurements made in thin and moderate cloud conditions (Toronto has about 44 % of daylight hours with clear-sky conditions per year). The relative strength of direct Sun measured by a collocated total sky imager (model TSI-880) is plotted at the top of each panel in Fig. 8 as an additional indicator of sky conditions. The relative strength of the direct Sun is from the integration of blocking-strip luminance. In general, when the relative strength of direct-Sun is high (> 60), good-quality direct-Sun and zenith-sky NO₂ data can both be produced. However, when Sun strength is moderate (30–60), only zenith-sky NO₂ data are reliable. When Sun strength is low (< 30), zenith-sky NO₂ has increased bias and needs to be filtered out.

5 Discussion

This study evaluated the performance of Pandora zenith-sky measurements with Pandora direct-Sun measurements, satellite measurements, and in situ measurements. In general, the quality of zenith-sky data is affected by three main factors: (1) quality of empirical zenith-sky AMFs; (2) cloud conditions (heavy clouds or moderate/thin clouds); and (3) quality of modelled NO₂ profile (this factor only applies to Pandora surface NO₂ data). The quality of empirical zenith-sky AMFs and the cloud effect have been addressed in Appendices A and C, respectively. The third factor is discussed in Sect. 5.1 and 5.2. The uncertainty estimations for Pandora zenith-sky and direct-Sun data products are provided in Appendix D.

5.1 Diurnal and seasonal variation

From the Pandora zenith-sky and direct-Sun measurements, and modelled NO₂ profiles, surface NO₂ concentrations were obtained that agree well with in situ measurements collected at the same location. The Pandora surface NO₂ data were also analyzed in more detail with a focus on temporal variations. Figure 9 shows the averaged surface NO₂ diurnal variations of four different datasets. The in situ instrument produces continuous measurements 24 h d⁻¹, whereas Pandora only has measurements when sunlight is available. The diurnal variation of surface NO₂ concentration is controlled by dynamics (e.g., vertical mixing, wind direction), photochemistry, and local emissions. Thus, the diurnal variations are calculated using only the hours when in situ, direct-Sun, and zenith-sky data are all available.

Figure 9 Diurnal variation of surface NO₂ concentration (2016–2017). The x axis is the local standard time (LST). Lines with dot/square symbols represent the hourly mean of corresponding data indicated by the legend. The shaded area represents the 1σ envelope.

Figure 9 shows that all four datasets/curves captured the enhanced morning surface NO₂ and the decreasing trend afterwards. However, the model has a positive offset (6–9 ppbv) in the morning (due in part to the use of older emission inventories; Moran et al., 2018) and a negative offset (1–3 ppbv) in the evening relative to the in situ data. For example, at 07:00 LST, in situ NO₂ is 14.9±9.3 ppbv, while GEM-MACH, Pandora DS, and Pandora ZS NO₂ are 23.5±15.0, 15.6±10.5, and 15.2±6.8 ppbv, respectively. At 17:00 LST, in situ NO₂ is 7.3±5.8 ppbv, while GEM-MACH, Pandora DS, and Pandora ZS NO₂ are 5.6±5.0, 3.6±2.6, and 5.2±3.4 ppbv, respectively. The larger standard deviations in the morning are due to the datasets not being divided into workdays and weekends. Compared to the modelled data, the Pandora direct-Sun and zenith-sky data show improvements in the morning but almost no changes for the evening. This feature is investigated and found to be correlated with the GEM-MACH modelled PBL height (details in Sect. 5.2).

Figure 10 Diurnal variation of surface NO₂ concentration by season (2016–2017). The x axis is the local standard time (LST). Each panel represents data collected in one season (spring, summer, autumn, or winter). Solid lines represent mean of corresponding data indicated by the legend. The shaded area represents the 1σ envelope.

The diurnal variation is also examined by grouping the data by seasons. Figure 10 shows that the surface NO₂ concentrations in winter (December, January, and February) are higher than the corresponding values in summer (June, July, and August). This difference is mainly due to short sunlit periods and less solar radiation (e.g., increased lifetime of NO₂ and decreased PBL height) in winter. The model has better agreement with the in situ data in summer than in the colder seasons. The best performance of the model is found around local noon, and this feature is not dependent on seasons. Figure 10 also shows that the quality of Pandora zenith-sky and direct-Sun surface NO₂ estimates is affected by the quality of GEM-MACH modelled data. For example, Fig. 10c shows that in autumn (September, October, and November), GEM-MACH has the largest offset in the morning. This error is thus propagated to the Pandora direct-Sun surface data, and leads to a larger offset in the morning (than any other season). On the other hand, when GEM-MACH shows a better agreement with in situ measurements (e.g., in spring and summer), Pandora zenith-sky and direct-Sun estimates also show better agreement with in situ observations. In general, both Pandora direct-Sun and zenith-sky surface NO₂ data show good agreement with in situ measurements in all seasons; the hourly mean values of Pandora surface NO₂ are all well within the 1σ envelope of the in situ measurements.

5.2 Planetary boundary layer effect

The larger morning offset in modelled surface NO₂ may indicate that the GEM-MACH modelled PBL heights are biased in the morning when the boundary layer is shallow. Figure 11 (left column) shows the modelled PBL height plotted as a function of the difference between modelled and in situ surface NO₂. Figure 11a shows that, in general, the difference between modelled and in situ NO₂ decreases with an increase of PBL height. When the modelled PBL height is less than 100 m, the mean difference is 18±12 ppbv (1σ), while when the modelled PBL height is 1 km, the mean difference is only 2.9±6.4 ppbv.

Figure 11 Illustration of planetary boundary layer (PBL) effect (2016–2017). The y axis is planetary boundary layer height in kilometers. The x axes for the left column are the difference between GEM-MACH and in situ surface NO₂ concentrations; the x axes for the right column are the difference between Pandora zenith-sky (C_pan-model) and in situ surface NO₂ concentration. Panels (a) and (b) show all available data, panels (c) and (d) show the morning data (before 09:00 LST), panels (e) and (f) show the noon data (from 11:00 to 13:59 LST), and panels (g) and (h) show the evening data (after 15:00 LST).

Even though the modelled surface concentrations are significantly impacted by the PBL, the modelled conversion ratio (from column to surface concentrations) seems unaffected since the surface NO₂ concentrations derived from Pandora zenith-sky data (C_pan-model) show much less dependence on the PBL height (Fig. 11b). When the modelled PBL height is less than 100 m, the mean difference is 0.9±8.9 ppbv. When the modelled PBL height is 1 km, the mean difference is slightly improved to 0.1±4.4 ppbv. Figure 11c and h show similar plots to Fig. 11a and b, but the dataset has been divided into three time bins (before 09:00, 11:00 to 13:59, and after 15:00 LST). Figure 11c, e, and f confirm that whenever the modelled PBL height is low, the relative difference between the model and in situ data is high. However, in general, most of these shallow PBL height conditions occur in the morning, and thus the modelled surface NO₂ has larger bias compared to in situ data in the morning. Figure 11d, f, and h show that Pandora zenith-sky surface NO₂ data have similar performance for all these three time bins, which indicates that the data have less PBL height dependency than the modelled data. In other words, the model is able to capture the ratio between the boundary layer partial column and surface NO₂, although the PBL height may not be correct in the model. When this ratio is applied to both Pandora direct-Sun and zenith-sky data, the estimated surface concentrations agree better with the in situ measurements.

6 Conclusions

The Pandora spectrometer was originally designed to retrieve total columns of trace gases such as ozone and NO₂ from direct-Sun spectral measurements in the UV–visible spectrum. In this work, a new zenith-sky total column NO₂ retrieval algorithm has been developed. The algorithm is based on empirical AMFs derived from nearly simultaneous direct-Sun and zenith-sky measurements. It is demonstrated that this algorithm can retrieve total columns in thin and moderate cloud conditions when direct-Sun measurements are not available: only 10 % of the measurements affected by heavy cloud have to be filtered out due to large systematic biases (68 %). The new Pandora zenith-sky total column NO₂ data shows only −4 % bias compared to the standard Pandora direct-Sun data product. In addition, OMI NO₂ SPv3.0 data demonstrate similar biases (−30 % and −29 %, respectively) when compared to direct-Sun and zenith-sky Pandora total column NO₂ data.

Surface NO₂ concentrations were calculated from Pandora direct-Sun and zenith-sky total column NO₂ using column-to-surface ratios derived from GEM-MACH regional chemical transport model. The bias between Pandora-based direct-Sun and zenith-sky NO₂ surface concentration estimates and in situ measurements is only −8 % and −7 % (with correlation coefficients 0.80 and 0.77), respectively, while the bias between the modelled concentrations and in situ measurements is up to 37 %. The Pandora-based surface NO₂ concentrations also show good diurnal and seasonal variation when compared to the in situ data. High surface NO₂ concentrations in the morning (from 06:00 to 09:00 LST) are present in all measured and modelled datasets, while, on average, the model overestimates surface NO₂ in the morning by 8.6 ppbv (at 07:00 LST). It appears that the bias in modelled surface NO₂ is related at least in part to an incorrectly diagnosed PBL height. In contrast, the difference between Pandora-based and in situ NO₂ does not show any significant dependence on the PBL height. Thus, to enable a fast and practical Pandora surface NO₂ data production, the use of a pre-calculated conversion ratio PSC-LUT is recommended.

The new retrieval algorithm for Pandora zenith-sky NO₂ measurements can provide high-quality NO₂ data (both total column and surface concentration) not only in clear-sky conditions but also in thin and moderate cloud conditions, when direct-Sun observations are not available. Long-term Pandora zenith-sky NO₂ data could be used in future satellite validation for the medium cloudy scenes. Moreover, a column-to-surface conversion look-up table was produced for the Pandora instruments deployed in Toronto; therefore, quick and practical Pandora-based surface NO₂ concentration data can be obtained for air quality monitoring purposes. The variation of conversion ratios in the PSC-LUT demonstrates that the surface NO₂ concentration is controlled not only by the PBL height but also by both boundary layer dynamics and photochemistry. This conversion approach can also be used to derive surface concentrations from satellite VCD measurements and thus can be particularly useful for the new generation of geostationary satellite instruments for air quality monitoring such as the Tropospheric Emissions: Monitoring of Pollution (TEMPO; Zoogman et al., 2014). Currently, the standard Pandora observation schedule includes direct-Sun, zenith-sky, and multi-axis scanning measurements (i.e., measuring at multiple viewing angles). At present, multi-axis measurement algorithms are still under development, but in the future, by using the multi-axis measurements and optimal estimation techniques (e.g., Rodgers, 2000) or the five-angle O₂O₂ ratio algorithm (Cede, 2019), it may be possible for Pandora measurements to be used to derive NO₂ tropospheric profiles and columns.

Appendix A: Empirical zenith-sky AMF

Before calculating the empirical zenith-sky AMF, the VCD_DS and dSCD_ZS have both been strictly filtered to ensure any measurements used in this calculation have the highest quality. For VCD_DS, data are filtered following Cede (2019) with several factors being considered, such as wavelength shift and residual in spectra fitting, direct-Sun AMF, and estimated uncertainties for the vertical column. For dSCD_ZS, data are filtered using similar criteria as for VCD_DS, with adjustments for zenith-sky observations.

The VCD_DS and dSCD_ZS data are merged and divided into several SZA bins. Each bin covers 5^∘. A multi-non-linear regression is performed by using the following equation:

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}\left[\begin{array}{c}{\mathrm{VCD}}_{\mathrm{1}}\\ \mathrm{⋮}\\ {\mathrm{VCD}}_n\end{array}\right]& =\left[\begin{array}{ccc}{\mathrm{dSCD}}_{\mathrm{1}}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\cdots }& {\mathrm{dSCD}}_n\end{array}\right]\left[\begin{array}{c}{b}_{\mathrm{1}}\\ \mathrm{⋮}\\ b_n\end{array}\right]\\ & +\mathrm{RCD}\left[\begin{array}{ccc}{I}_{\mathrm{1}}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\cdots }& I_n\end{array}\right]\left[\begin{array}{c}{b}_{\mathrm{1}}\\ \mathrm{⋮}\\ b_n\end{array}\right],\end{array}\end{array}$$

where VCD_n is not a single direct-Sun VCD data point but is an m×1 matrix (m is the total number of measurements in SZA bin number n); the VCD_n represents all direct-Sun VCDs in a 5^∘ SZA bin, and each element of the m×1 matrix is a single VCD in that SZA bin. Similarly, dSCD_n is also an m×1 matrix, with each element representing a single coincident zenith-sky dSCD in SZA bin number n. I_n is an m×1 indicator function, where the elements of I_n are set to 1. The RCD and b₁ to b_n are the parameters to be retrieved. In short, the design of this regression is based on Eq. (2) (Sect. 3.1). The idea is to retrieve zenith-sky AMFs in several SZA bins, and, at the same time, all these regressions in different SZA bins are constrained to share a common predictor (RCD). The regression model can be solved by using an iterative procedure (Seber and Wild, 2003) to yield the estimated coefficients, b₁ to b_n and RCD. The b_n is the reciprocal of zenith-sky AMF in SZA bin n.

This regression model has been evaluated by using different sizes for the SZA bins. A 5^∘ SZA bin is selected because the SZA bin must be small enough to capture the SZA dependency on zenith-sky AMFs, and, at the same time, it must also be large enough to ensure a sufficient number of measurements in each SZA bin (to perform reliable regressions). In order to deal with the diurnal variation of NO₂ concentration and changing of profile shape (e.g., due to changing of boundary layer heights), the dataset has been divided into morning and evening sets, and discrete AMFs are retrieved for morning and evening hours separately (see the blue and red squares with error bars in Fig. 1).

Next, these discrete AMF values are used to fit an empirical zenith-sky NO₂ AMF function, which has the expression

$$\begin{array}{}\text{(A2)}& \mathrm{AMF}=a_{\mathrm{1}}+\left(\mathrm{1.02}-a_{\mathrm{1}}\right)/\cos \left(\mathrm{SZA}\right).\end{array}$$

The fitted empirical zenith-sky AMFs are shown in Fig. 1 as blue and red lines (data regression period from February 2015 to September 2017). Several sensitivity tests have been performed to assess the quality of the empirical zenith-sky AMFs, including fitting the AMFs with/without a diurnal difference, fitting the AMFs with different empirical functions (e.g., exponential and simple geometry approximation), and fitting the AMFs by seasons. All these different choices of empirical AMFs fitting functions or methods only introduce less than 5 % difference in the retrieved empirical AMFs. Thus, to make the empirical AMFs simple and robust, we selected to fit with a diurnal difference (Eq. 5). In addition, the current empirical AMFs are limited to high- and intermediate-Sun conditions (i.e., SZA < 75^∘). For low-Sun conditions, the total AMF for zenith-sky measurements is expected to be a strong function of not only the SZA but also the tropospheric column itself. Thus, for future work to derive low-Sun empirical zenith-sky AMFs, the stronger influence of PBL NO₂ has to be accounted (i.e., the geometry-form AMFs are not enough).

Appendix B: Stratospheric NO₂ column

Several stratospheric NO₂ column values were tested and used in the surface NO₂ concentration algorithm (Eq. 3). Figure B1a shows the OMI monthly mean (referred to as OMI) and Pratmo box-model stratospheric column NO₂ (Adams et al., 2016; McLinden et al., 2000) (referred to as box). Since the satellite only samples Toronto once per day, the OMI stratospheric NO₂ lacks diurnal variation. To account for the diurnal variation, diurnal ratios of NO₂ VCD have been calculated and applied to OMI monthly mean data. The stratospheric NO₂ columns are calculated using

$$\begin{array}{}\text{(B1)}& V_{\mathrm{OMI}}\left(t\right)={\displaystyle \frac{V_{\mathrm{box}}\left(t\right)}{V_{\mathrm{box}}\left(t_{\mathrm{0}}\right)}}\times V_{\mathrm{OMI}}\left(t_{\mathrm{0}}\right),\end{array}$$

where V_OMI(t₀) is the OMI-measured stratospheric column, t₀ is OMI overpass time, V_box(t₀) is the modelled stratospheric column at OMI overpass time, V_box(t) is the modelled stratospheric column at time t, and V_OMI(t) is the interpolated stratospheric column at time t. The interpolated OMI stratospheric columns are referred to as OMI box. The grey dots in Fig. B1b are OMI-box stratospheric NO₂ columns. The monthly mean of the box model (blue line) and OMI box (black line) show that the amplitude of OMI box is larger than the amplitude of the box model.

Figure B1 Time series of measured and modelled NO₂ columns: (a) stratospheric columns from the box model (hourly) and OMI (monthly), (b) stratospheric columns from OMI box (hourly), box (monthly) and OMI box (monthly), and (c) total columns from Pandora zenith sky and OMI.

To justify why this diurnal variation has to be included, Fig. B1c shows the total column NO₂ time series. The diurnal stratospheric NO₂ variation is about 0.1 DU in the summer (see grey dots in Fig. A1b) when Pandora measured monthly mean total column is about 0.5 DU (Fig. B1c). Thus, neglecting this diurnal variation will lead to diurnal biases in the derived surface NO₂ data (e.g., in the morning, this will lead to the overestimation of the stratospheric NO₂ and thus the underestimation of surface NO₂). Please note that the strength of this bias is related to (1) the NO₂ profile (weights between stratospheric and tropospheric NO₂), and (2) the observation geometry (direct Sun or zenith sky). In general, an urban site with direct-Sun observation should have smaller impact from the stratospheric diurnal variation. On the other hand, a rural site with zenith-sky observation should have a significant impact.

Appendix C: Cloud effect and heavy cloud filtration

Direct-Sun measurements need an unobscured Sun. Even thin clouds could decrease the quality of retrieved NO₂ total columns, especially for low-altitude clouds. Unlike direct-Sun measurements, zenith-sky observations are made with scattered sunlight and have limited sensitivity to cloud cover. For example, Hendrick et al. (2011) calculated that, for NDACC UV–visible zenith-sky ozone measurements, clouds only contribute 3.3 % to the total random error. This is because a trace gas that is mostly distributed in the stratosphere has the mean scattering layer located at a higher altitude than the cloud layer. However, this assumption may not be valid for NO₂. Depending on the properties of the clouds and the NO₂ profile, the clouds could have non-negligible impacts on zenith-sky NO₂ observations.

A typical method of removing zenith-sky measurements affected by heavy clouds is to eliminate measurements with large enhancements of O₄ and/or H₂O (Van Roozendael and Hendrick, 2012). In the Pandora zenith-sky NO₂ retrieval, we use the O₄ dSCDs. Since the measured O₄ dSCDs has SZA dependency, all measured O₄ dSCDs are plotted against SZA, and a second-order quantile regression (Koenker and Hallock, 2001) is applied to select the top few percentiles of the measured O₄ dSCDs.

Figure C1 Illustration of cloud effect and heavy cloud data filtration: panel (a) shows measured O₄ differential slant column densities vs. solar zenith angle; the grey dots represent the top 0–10th percentile range of O₄. Panel (b) shows the scatter plot of zenith-sky vs. in situ surface NO₂ data that has O₄ values within the 0–10th percentile range (as identified in panel a); panel (c) is similar to (a) but the grey dots represent the 40th–50th percentile range of O₄; panel (d) is similar to (b) but uses the data that has O₄ value within the 40th–50th percentile range. On the scatter plots in panels (b) and (d), the red line is the linear fit with intercept set to 0 and the black line is the one-to-one line. All plots are colour-coded by the normalized density of the points.

Figure C2 Correlation coefficient and bias (slope) between zenith-sky and in situ surface NO₂ data in different O₄ dSCD percentile bins. Panel (a) shows the correlation coefficients; panel (b) shows the slopes of linear fit with intercept set to 0.

Figure C3 Example of surface NO₂ concentration time series in all conditions. The in situ, Pandora DS, and Pandora ZS surface NO₂ concentrations are shown by different coloured dots. The TSI relative strength of direct-Sun data is plotted as a colour-coded horizontal dotted line in the top area of each panel. For Pandora zenith-sky data, the measurements with enhanced O₄ (heavy cloud indicator) are also labelled by green squares. Dates are in mm/dd format.

Figure C1 shows examples of selected Pandora zenith-sky NO₂ data and their corresponding O₄ dSCD values. For example, Fig. C1a shows the O₄ dSCDs vs. SZA, and the top 10 percentiles of the data with enhanced O₄ are marked in grey. The corresponding Pandora zenith-sky data are plotted against in situ data in Fig. C1b, which shows low correlation (R=0.34) and high bias (68 %). This result indicates that the enhanced scattering due to heavy clouds caused a positive bias in the Pandora zenith-sky NO₂ retrieval. Figures C1c and d are similar to Fig. C1a and b but for selected Pandora zenith-sky NO₂ data that have O₄ values within the 40th to 50th percentile range. Figure C1d shows that when O₄ is not enhanced, the derived zenith-sky NO₂ has good agreement with in situ data (R=0.8 and bias $=-\mathrm{5}$ %). To summarize how the retrieved O₄ dSCDs can indicate the quality of the Pandora zenith-sky NO₂, the results from the other percentile bins are shown in Fig. C2. In general, besides the top 10th percentile of data, the results from all the other bins show good correlation (above 0.6) and low bias. Thus, in this study, the Pandora zenith-sky NO₂ data that have O₄ values within only the top 10th percentile are considered to be affected by heavy clouds and are removed. Some examples of this heavy cloud effect are shown in Figs. C3 and 8 in Sect. 4.3.

Appendix D: Uncertainty estimation

The uncertainties of retrieved Pandora zenith-sky NO₂ data products (total column and surface concentration) are estimated and discussed here to assess the quality of the data products. The uncertainties of total column and surface concentrations are estimated first using the uncertainty propagation method (referred to here as the UP method) based on Eqs. (2) and (3). The combined uncertainties of total column can be calculated using

$$\begin{array}{}\text{(D1)}& {\mathit{\sigma }}_{{\mathrm{VCD}}_{\mathrm{ZS}}}=\sqrt[\mathrm{2}]{{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{dSCD}}}{\mathrm{AMF}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{RCD}}}{\mathrm{AMF}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{AMF}}\times \mathrm{SCD}}{{\mathrm{AMF}}^{\mathrm{2}}}}\right)}^{\mathrm{2}}},\end{array}$$

where σ_dSCD is the statistical uncertainty on the DOAS fit (output of QDOAS) and σ_RCD and σ_AMF are the estimated statistical uncertainties using standard errors of the RCD and the zenith-sky empirical AMF regression, respectively (Eq. A1). To estimate the upper limit of the nominal uncertainty, AMF and SCD are used as median and maximum values in the dataset, respectively.

The combined uncertainties of the surface concentration can be calculated using

$$\begin{array}{}\text{(D2)}& \begin{array}{rl}{\mathit{\sigma }}_{C_{\mathrm{Pan}}}& =\sqrt[\mathrm{2}]{{\left(R_{\mathrm{CV}}{\mathit{\sigma }}_{V_{\mathrm{Pan}}}\right)}^{\mathrm{2}}+{\left(R_{\mathrm{CV}}{\mathit{\sigma }}_{V_{\mathrm{strat}}}\right)}^{\mathrm{2}}}\\ & \stackrel{\mathrm{‾}}{+{\left(R_{\mathrm{CV}}{\mathit{\sigma }}_{V_{\mathrm{ftrop}}}\right)}^{\mathrm{2}}+{\left(V_{\mathrm{Pan}}-V_{\mathrm{strat}}-V_{\mathrm{ftrop}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{R}}^{\mathrm{2}}},\end{array}\end{array}$$

where ${\mathit{\sigma }}_{V_{\mathrm{pan}}}$ is the uncertainty of Pandora zenith-sky total column NO₂, (here we use the derived σ_VCD in Eq. D1), ${\mathit{\sigma }}_{V_{\mathrm{strat}}}$ is the uncertainty of the stratospheric NO₂ column (estimated using the 1σ standard deviation of the V_strat), ${\mathit{\sigma }}_{V_{\mathrm{ftrop}}}$ is the uncertainty of the free troposphere NO₂ column (estimated using the 1σ standard deviation of the V_ftrop). R_CV is the GEM-MACH calculated surface VMR to PBL column ratio, and σ_R is the uncertainty of that ratio (estimated using the 1σ standard deviation of the R_CV). The means of R_CV, V_Pan, V_strat, and V_ftrop are used in the uncertainty estimation.

Table D1 Estimated uncertainties for Pandora zenith-sky total column and surface NO₂.

Besides the UP method, another simple approach to estimate uncertainty is to compare the data product with another high-quality (lower uncertainty) coincident datum. For example, if we assume that the Pandora direct-Sun total column NO₂ data can represent the true value, we can estimate the uncertainty of Pandora zenith-sky total column NO₂ by calculating the 1σ standard deviation of their difference (referred to here as the SDD method):

$$\begin{array}{}\text{(D3)}& {\mathit{\sigma }}_{{\mathrm{VCD}}_{\mathrm{ZS}}}=\mathit{\sigma }\left({\mathrm{VCD}}_{\mathrm{DS}}-{\mathrm{VCD}}_{\mathrm{ZS}}\right).\end{array}$$

Similarly, if we assume that the in situ surface NO₂ VMR can represent the true value, the uncertainty of Pandora zenith-sky-based surface NO₂ VMR can be given by

$$\begin{array}{}\text{(D4)}& {\mathit{\sigma }}_{C_{\mathrm{Pan}}}=\mathit{\sigma }\left(C_{\text{in situ}}-C_{\mathrm{pan}}\right).\end{array}$$

Also, if there is systematic bias between the two datasets, it can be removed and the random uncertainty can be calculated by

$$\begin{array}{}\text{(D5)}& {\displaystyle }{\mathit{\sigma }}_{{\mathrm{VCD}}_{\mathrm{ZS}}}=\mathit{\sigma }\left({\mathrm{VCD}}_{\mathrm{DS}}-k_{\mathrm{1}}{\mathrm{VCD}}_{\mathrm{ZS}}\right),\text{(D6)}& {\displaystyle }{\mathit{\sigma }}_{C_{\mathrm{Pan}}}=\mathit{\sigma }\left(C_{\text{in situ}}-k_{\mathrm{2}}{C}_{\mathrm{pan}}\right),\end{array}$$

where k₁ and k₂ are the slopes in the linear fits with intercept set to zero (e.g., slopes in Figs. 2 and 6). This method is referred to here as the unbiased SDD. These three uncertainty estimation methods (UP, SDD, and unbiased SDD) were all implemented, and the results are summarized in Table D1. The results show that Pandora zenith-sky total column NO₂ data have a 0.09–0.12 DU uncertainty that is about twice the Pandora direct-Sun total column nominal accuracy (0.05 DU, at 1σ level). When using the UP method, for the worst-case scenario, the Pandora zenith-sky total column NO₂ data have a 0.17 DU uncertainty (i.e., using minimum of AMFs to estimate the upper limit of uncertainty). The estimated Pandora zenith-sky-based surface NO₂ VMR data have uncertainties from 4.8 to 6.5 ppbv. In Eq. (D2), the contributions of the V_Pan, V_Strat, V_ftrop, and R_CV terms to the total uncertainty are 36 %, 2 %, 0.3 %, and 62 %, respectively. This result indicates that the uncertainty in the Pandora zenith-sky-based surface NO₂ VMR is dominated by the uncertainties of Pandora zenith-sky total column NO₂ and the modelled column-to-surface conversion ratio (R_CV). However, note that this uncertainty budget depends on the NO₂ vertical distributions and hence may vary from site to site; e.g., in Toronto, tropospheric column NO₂ is typically 2–4 times higher than stratospheric column NO₂, and thus the contribution to uncertainty from V_Pan is much larger than the corresponding contributions from V_Strat and V_ftrop. In addition, the uncertainty of Pandora direct-Sun surface NO₂ VMR is also estimated and provided in Table D1. It shows slightly better results than for zenith-sky-based surface NO₂ VMR.