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

Abstract

Pandora spectrometers can retrieve nitrogen dioxide (${\mathrm{NO}}_{\mathrm{2}}$) vertical column densities (VCDs) via two viewing geometries: direct Sun and zenith sky. The direct-Sun ${\mathrm{NO}}_{\mathrm{2}}$ 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${\mathrm{NO}}_{\mathrm{2}}$ measurements in cloudy conditions, a simple but robust ${\mathrm{NO}}_{\mathrm{2}}$ retrieval algorithm is developed for Pandora zenith-sky measurements. This algorithm derives empirical zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ AMFs from coincident high-quality direct-Sun ${\mathrm{NO}}_{\mathrm{2}}$ observations. Moreover, the retrieved Pandora zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ VCD data are converted to surface ${\mathrm{NO}}_{\mathrm{2}}$ concentrations with a scaling algorithm that uses chemical-transport-model predictions and satellite measurements as inputs. ${\mathrm{NO}}_{\mathrm{2}}$ VCDs and surface concentrations are retrieved from Pandora zenith-sky measurements made in Toronto, Canada, from 2015 to 2017. The retrieved Pandora zenith-sky${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ data also agree well with in situ measurements (diurnal difference within $\pm \mathrm{2}$ ppbv). Overall, this work shows that the new Pandora zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ 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 (${\mathrm{NO}}_{\mathrm{2}}$) 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. ${\mathrm{NO}}_{\mathrm{2}}$ is a nitrate aerosol precursor, and it also contributes to acid deposition and eutrophication (ECCC, 2016). Exposure to ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ concentrations are regulated by many environmental agencies (e.g., Environment and Climate Change Canada and US Environment Protection Agency), in situ ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ air quality monitoring is the photolytic conversion/chemiluminescence technique, which converts ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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${}^{-\mathrm{2}}$ or Dobson units (1 DU $=\mathrm{2.6870}\times {\mathrm{10}}^{\mathrm{16}}$ molec cm${}^{-\mathrm{2}}$). The Pandora direct-Sun ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ VCD data products available. As one goal of this work, we have focused on developing a new ${\mathrm{NO}}_{\mathrm{2}}$ retrieval algorithm for zenith-sky measurements to expand Pandora ${\mathrm{NO}}_{\mathrm{2}}$ measurements into cloudy scenes.

In addition to retrieval of zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$, another goal of this work is to derive surface ${\mathrm{NO}}_{\mathrm{2}}$ concentration from total column measurements. Surface ${\mathrm{NO}}_{\mathrm{2}}$ has been a focus of scientific studies due to its strong correlation with air quality (AQ) and health issues (ECCC, 2016), with ${\mathrm{NO}}_{\mathrm{2}}$ as one of the three components (along with ozone and PM${}_{\mathrm{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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ surface mixing ratios from Pandora direct-Sun total column ${\mathrm{NO}}_{\mathrm{2}}$ via application of a planetary boundary layer (PBL) height correction factor. Kollonige et al. (2017) adapted this method and compared Pandora direct-Sun surface ${\mathrm{NO}}_{\mathrm{2}}$ and Ozone Monitoring Instrument (OMI) surface ${\mathrm{NO}}_{\mathrm{2}}$. They concluded that the two main sources of error for the conversion of the total column ${\mathrm{NO}}_{\mathrm{2}}$ to surface ${\mathrm{NO}}_{\mathrm{2}}$ are (1) poor weather conditions (e.g., cloud cover and precipitation) and (2) PBL height estimation, both of which affect the ${\mathrm{NO}}_{\mathrm{2}}$ column–surface relationship and instrument sensitivities to boundary layer ${\mathrm{NO}}_{\mathrm{2}}$. Thus, in this work, we present a simple but robust algorithm for deriving surface ${\mathrm{NO}}_{\mathrm{2}}$ concentration from Pandora zenith-sky measurements, which has several advantages, such as the ability (1) to extend Pandora ${\mathrm{NO}}_{\mathrm{2}}$ measurements to cloudy conditions and (2) to provide more accurate surface ${\mathrm{NO}}_{\mathrm{2}}$ concentration estimates that are less sensitive to PBL height. This work also provides reliable total column ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ data used in this study. In Sect. 3, the empirical AMFs for Pandora zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ measurements are derived using high-quality Pandora direct-Sun total column ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ retrievals to help evaluate the performance of the empirical AMFs. Also, the retrieved Pandora zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$ data are evaluated by comparison with satellite measurements. In Sect. 4, the zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$ data are converted to surface concentration by using a scaling algorithm. The zenith-sky surface ${\mathrm{NO}}_{\mathrm{2}}$ concentration data are assessed by comparison with in situ measurements. Lastly, in Sect. 5, several aspects of this zenith-sky surface ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$

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 $\mathrm{\mu }$m entrance slit, 1200 groove mm${}^{-\mathrm{1}}$ grating, and a $\mathrm{2048}\times \mathrm{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, ${\mathrm{NO}}_{\mathrm{2}}$, and sulfur dioxide (${\mathrm{SO}}_{\mathrm{2}}$), are fitted to the spectra.

The Pandora direct-Sun total column ${\mathrm{NO}}_{\mathrm{2}}$ data are produced using Pandora's standard ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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, $\mathrm{2}\mathit{\sigma }$ 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${}^{\circ }$ N, $-\mathrm{79.4680}$${}^{\circ }$ 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 ${\mathrm{NO}}_{\mathrm{2}}$

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 ${\mathrm{NO}}_{\mathrm{2}}$. Following the NDACC recommendations (Van Roozendael and Hendrick, 2012), ${\mathrm{NO}}_{\mathrm{2}}$ differential slant column densities (dSCDs) are retrieved in the 425–490 nm window (to retrieve oxygen collision complex simultaneously). The oxygen collision complex (${\mathrm{O}}_{\mathrm{2}}$)${}_{\mathrm{2}}$ (referred here as ${\mathrm{O}}_{\mathrm{4}}$), 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). ${\mathrm{O}}_{\mathrm{4}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ at an effective temperature of 254.5 K (Vandaele et al., 1998), ozone at an effective temperature of 223 K (Bogumil et al., 2003), ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ (Rothman et al., 2005), ${\mathrm{O}}_{\mathrm{4}}$ (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 ${\mathrm{NO}}_{\mathrm{2}}$ dSCDs, which can be converted to total column ${\mathrm{NO}}_{\mathrm{2}}$ via the Langley plot method with the use of the NDACC ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$. Other empirical zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ AMFs have been developed and are used to convert ${\mathrm{NO}}_{\mathrm{2}}$ 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 $\mathrm{780}\times \mathrm{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 $\mathrm{13}\mathrm{km}\times \mathrm{24}\mathrm{km}$ (along track $\times$ across track), whereas those close to the track edge (e.g., view zenith angle $=\mathrm{56}$${}^{\circ }$) have a ground footprint roughly of $\mathrm{23}\mathrm{km}\times \mathrm{126}\mathrm{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 $\mathrm{13}\mathrm{km}\times \mathrm{24}\mathrm{km}$.

The OMI ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 $\times$ longitude) resolution a priori ${\mathrm{NO}}_{\mathrm{2}}$ 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, ${\mathrm{NO}}_{\mathrm{2}}$, ${\mathrm{SO}}_{\mathrm{2}}$, 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ concentration is measured using a photolytic ${\mathrm{NO}}_{\mathrm{2}}$ instrument (Thermo 42i) that is also sensitive to other gaseous inorganic nitrogen compounds (e.g., nitric acid (${\mathrm{HNO}}_{\mathrm{3}}$) and peroxyacetyl nitrate (PAN)) (McLinden et al., 2014). Thus, in areas where direct ${\mathrm{NO}}_x$ (nitrogen oxides) emission sources are limited and other nitrogen compounds are present, ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ concentrations (< 5 ppbv) (Yushan Su, Ontario Ministry of the Environment, Conservation and Parks, personal communication, October 2018).

2.2 Numerical models

Predicted ${\mathrm{NO}}_{\mathrm{2}}$ fields from three atmospheric chemistry models are used in the algorithm described in Sect. 4.1 to derive surface ${\mathrm{NO}}_{\mathrm{2}}$ concentration from Pandora zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ and surface ${\mathrm{NO}}_{\mathrm{2}}$ concentration. The stratospheric ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_x$ emissions in the PBL are included in the operational model; free-tropospheric ${\mathrm{NO}}_x$ emissions from lightning and in-flight aircraft are not considered. In this work, the GEM-MACH hourly ${\mathrm{NO}}_{\mathrm{2}}$ 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 $\mathrm{10}\mathrm{km}\times \mathrm{10}\mathrm{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 ${\mathrm{NO}}_{\mathrm{2}}$ (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 ${\mathrm{NO}}_x$ source of 6 Tg N yr${}^{-\mathrm{1}}$ (Martin et al., 2002) was imposed. Lightning ${\mathrm{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 $\times$ longitude) grid in nested mode over North America and was driven by assimilated meteorology from the Goddard Earth Observing System (GEOS-5). The modelled ${\mathrm{NO}}_{\mathrm{2}}$ profiles were used to calculate monthly mean ${\mathrm{NO}}_{\mathrm{2}}$ partial columns in the free troposphere (1.5 to 12 km), as the GEM-MACH model does not include free-tropospheric ${\mathrm{NO}}_{\mathrm{2}}$ 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 (${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$), methane (${\mathrm{CH}}_{\mathrm{4}}$), and water vapor (${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$)) and halogen families (${\mathrm{NO}}_y$, ${\mathrm{Cl}}_y$, and ${\mathrm{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 ${\mathrm{NO}}_{\mathrm{2}}$ has a strong diurnal variation; therefore, diurnal corrections must be applied when OMI stratospheric ${\mathrm{NO}}_{\mathrm{2}}$ measurements (around local noon) are interpolated to Pandora measurement times. Ratios of modelled stratospheric ${\mathrm{NO}}_{\mathrm{2}}$ columns are calculated at OMI overpass time and Pandora measurement time. These ratios are multiplied by the OMI measured stratospheric ${\mathrm{NO}}_{\mathrm{2}}$ to produce stratospheric ${\mathrm{NO}}_{\mathrm{2}}$ columns corresponding to the time of Pandora measurements. Details about the use of the Pratmo box model and the calculation of stratospheric ${\mathrm{NO}}_{\mathrm{2}}$ partial columns are provided in Sect. 4.1.

3

Total column ${\mathrm{NO}}_{\mathrm{2}}$ retrieval

3.1 Zenith-sky air mass factor

The NDACC UV–visible network uses zenith-sky AMFs in its total column ${\mathrm{NO}}_{\mathrm{2}}$ retrievals. To improve the overall homogeneity of the UV–visible ${\mathrm{NO}}_{\mathrm{2}}$ column measurements, NDACC recommended using the ${\mathrm{NO}}_{\mathrm{2}}$ AMF LUT (Van Roozendael and Hendrick, 2012). This LUT is based on climatological ${\mathrm{NO}}_{\mathrm{2}}$ profiles that are composed of (1) 20–60 km ${\mathrm{NO}}_{\mathrm{2}}$ profiles developed by Lambert et al. (1999, 2000) and (2) 12–20 km ${\mathrm{NO}}_{\mathrm{2}}$ profiles derived from SAOZ (Système D'Analyse par Observations Zénithales) balloon observations (Van Roozendael and Hendrick, 2012). The ${\mathrm{NO}}_{\mathrm{2}}$ concentration is set to zero below 12 km altitude. The ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ retrievals. Note that the absence of tropospheric ${\mathrm{NO}}_{\mathrm{2}}$ in the NDACC LUT construction will lead to an underestimation of the total column ${\mathrm{NO}}_{\mathrm{2}}$ in urban areas. For example, from 2015 to 2017, tropospheric ${\mathrm{NO}}_{\mathrm{2}}$ accounted for $\mathrm{73}\pm \mathrm{11}$ % ($\mathrm{1}\mathit{\sigma }$) of the total column amounts in Toronto (OMI SPv3.0 data). To account for this significant tropospheric ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ measurements, this work took a different but simple and robust approach to derive zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$.

Empirical AMFs are calculated for Pandora zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ measurements in such a way that they can be used to retrieve zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$ values that match the high-quality Pandora direct-Sun total column ${\mathrm{NO}}_{\mathrm{2}}$ 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 ($\mathit{\mu }$) 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 ${\mathrm{NO}}_{\mathrm{2}}$ dSCDs.

The relation between VCD and dSCD can be expressed as

1$$\mathrm{VCD}={\displaystyle \frac{\mathrm{dSCD}+\mathrm{RCD}}{\mathrm{AMF}}},$$ 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${}_{\text{ZS-Emp}}$) can be calculated by assuming VCD${}_{\mathrm{DS}}=$ VCD${}_{\mathrm{ZS}}$, which gives 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)}}.$$ Next, we can use nearly coincident VCD${}_{\mathrm{DS}}$ and dSCD${}_{\mathrm{ZS}}$ in a multi-non-linear regression to retrieve AMF${}_{\text{ZS-Emp}}$ and RCD${}_{\mathrm{ZS}}$ together. To ensure the quality of the retrieved AMF${}_{\text{ZS-Emp}}$, only high-quality direct-Sun total column ${\mathrm{NO}}_{\mathrm{2}}$ data are used with SZA < 75${}^{\circ }$. 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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${}_{\text{ZS-Emp}}$ and the one based on NDACC AMFs is referred to as VCD${}_{\text{ZS-NDACC}}$. The RCD value used in the retrievals is $\mathrm{0.39}\pm \mathrm{0.01}$ DU, which is retrieved along with AMF${}_{\text{ZS-Emp}}$ (Appendix A). Figure 2 shows the comparisons of the ${\mathrm{NO}}_{\mathrm{2}}$ 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 $\pm \mathrm{5}$ min of Pandora zenith-sky measurement, (2) SZA < 75${}^{\circ }$, and (3) Pandora direct-Sun total column ${\mathrm{NO}}_{\mathrm{2}}$ data have assured high quality (BlickP L2 data quality flag for nitrogen dioxide is 0). In general, the VCD${}_{\text{ZS-Emp}}$ and VCD${}_{\text{ZS-NDACC}}$ performed as expected. Compared with VCD${}_{\mathrm{DS}}$, the VCD${}_{\text{ZS-NDACC}}$ shows a $-\mathrm{25}$ % bias, while the VCD${}_{\text{ZS-Emp}}$ only shows a $-\mathrm{4}$ % bias (indicated by the red lines on each panel and their slopes). In addition, VCD${}_{\text{ZS-Emp}}$ shows less SZA dependence than VCD${}_{\text{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 ${\mathrm{NO}}_{\mathrm{2}}$ profile should be included when calculating empirical zenith-sky AMFs. In the rest of the paper, only the zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ retrieved using empirical AMFs will be discussed. The derived zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ data discussed in Sect. 3 are a “clear-sky subset” of Pandora zenith-sky measurements. The assessment of Pandora zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ measurements in cloudy conditions is provided in Sect. 4.

Figure 1

Comparison of zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ air mass factors. Blue and red squares with error bars (standard error) represent the empirical discrete zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ AMFs. NDACC AMFs calculated using the NDACC look-up table and assuming no ${\mathrm{NO}}_{\mathrm{2}}$ in the troposphere are shown in yellow.

[Figure omitted. See PDF]

Figure 2

Comparisons of ${\mathrm{NO}}_{\mathrm{2}}$ total columns (2015–2017): (a) zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$ retrieved using empirical AMFs vs. direct-Sun total column ${\mathrm{NO}}_{\mathrm{2}}$, (b) zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$ retrieved using NDACC AMFs vs. direct-Sun total column ${\mathrm{NO}}_{\mathrm{2}}$. 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).

[Figure omitted. See PDF]

3.2 Comparison with satellite measurements

To illustrate the ${\mathrm{NO}}_{\mathrm{2}}$ variability over Toronto, Fig. 3 shows the time series (2015–2017) from Pandora direct-Sun, zenith-sky, and OMI SPv3.0 total column ${\mathrm{NO}}_{\mathrm{2}}$. In general, the ${\mathrm{NO}}_{\mathrm{2}}$ datasets from the ground-based Pandora instrument and the satellite follow the same pattern. However, the satellite data are likely to miss the peak ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ in Toronto from 2015 to 2017.

[Figure omitted. See PDF]

We also performed regression analyses by using the following coincidence criteria: (1) nearest (in time) measurement that was within $\pm \mathrm{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 < $=\mathrm{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 $=\mathrm{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 ${\mathrm{NO}}_{\mathrm{2}}$, respectively. Figure 4c and d show similar comparisons but only use OMI ${\mathrm{NO}}_{\mathrm{2}}$ measured by “small pixels” (i.e., having viewing zenith angle of less than 35${}^{\circ }$). 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 ${\mathrm{NO}}_{\mathrm{2}}$, both OMI and Pandora zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$ underestimate the local ${\mathrm{NO}}_{\mathrm{2}}$ 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 $\mathrm{0.64}\pm \mathrm{0.02}$, $\mathrm{0.67}\pm \mathrm{0.02}$, $\mathrm{0.70}\pm \mathrm{0.04}$, and $\mathrm{0.71}\pm \mathrm{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 ${\mathrm{NO}}_{\mathrm{2}}$, and the relatively large OMI pixels are also generally sampling areas of less ${\mathrm{NO}}_{\mathrm{2}}$ in the vicinity. The use of the relatively coarse (1${}^{\circ }$) GMI model for profile shapes (Sect. 2.1.3) will also lead to a low bias considering the peak ${\mathrm{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 ${\mathrm{NO}}_{\mathrm{2}}$ (2015–2017). Panels (a) and (c) show OMI vs. Pandora direct-Sun ${\mathrm{NO}}_{\mathrm{2}}$, and (b) and (d) show OMI vs. Pandora zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$. 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.

[Figure omitted. See PDF]

Ialongo et al. (2016) reported a similar negative bias using OMI SPv3.0 and Pandora direct-Sun total column ${\mathrm{NO}}_{\mathrm{2}}$ in Helsinki ($-\mathrm{32}$ % bias and $R=\mathrm{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 ${\mathrm{NO}}_{\mathrm{2}}$ 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=\mathrm{0.3}$). The authors attributed the poor agreement to the coarse resolution of OMI and its AMFs computed with GMI a priori ${\mathrm{NO}}_{\mathrm{2}}$ profiles. In general, our comparison results show that (1) the Pandora direct-Sun total column ${\mathrm{NO}}_{\mathrm{2}}$ data measured in Toronto have a reasonable agreement with OMI, and (2) the Pandora zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$ data show results similar to those for direct-Sun total column when compared with OMI SPv3.0.

4

Surface ${\mathrm{NO}}_{\mathrm{2}}$ concentration retrieval

The performance of the clear-sky Pandora zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ concentration from Pandora zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$ measurements. Following Lamsal et al. (2008) and McLinden et al. (2014), the surface ${\mathrm{NO}}_{\mathrm{2}}$ concentration is estimated using the modelled profile and surface concentration:

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}},$$ where $C_{\mathrm{pan}}$ is the surface ${\mathrm{NO}}_{\mathrm{2}}$ volume mixing ratio (VMR) to be estimated, $C$ is the surface ${\mathrm{NO}}_{\mathrm{2}}$ VMR from GEM-MACH (or G-M), $V_{\mathrm{pan}}$ is the total column ${\mathrm{NO}}_{\mathrm{2}}$ measured by Pandora, $V_{\mathrm{strat}}$ is the stratospheric ${\mathrm{NO}}_{\mathrm{2}}$ partial column, $V_{\mathrm{ftrop}}$ is the ${\mathrm{NO}}_{\mathrm{2}}$ partial column in the free troposphere, and $V_{\mathrm{PBL}}$ is the ${\mathrm{NO}}_{\mathrm{2}}$ 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_{\mathrm{PBL}}$ (0–1.5 km) is integrated from the GEM-MACH ${\mathrm{NO}}_{\mathrm{2}}$ profile and $V_{\mathrm{ftrop}}$ (1.5–12 km) is integrated from the GEOS-Chem ${\mathrm{NO}}_{\mathrm{2}}$ profile. Both GEM-MACH and GEOS-Chem have an hourly temporal resolution. Thus, the integrated $V_{\mathrm{PBL}}$ and $V_{\mathrm{ftrop}}$ can account for ${\mathrm{NO}}_{\mathrm{2}}$ diurnal variation. However, $V_{\mathrm{strat}}$ is from OMI monthly mean stratospheric ${\mathrm{NO}}_{\mathrm{2}}$, which does not have diurnal variation. Thus, the Pratmo box model is used to calculate stratospheric ${\mathrm{NO}}_{\mathrm{2}}$ diurnal ratios. The OMI stratospheric ${\mathrm{NO}}_{\mathrm{2}}$ columns are interpolated to morning and evening hours by multiplying by the box-model diurnal ratios. Details about the calculation of $V_{\mathrm{strat}}$ as well as references are provided in Appendix B.

The $(C/V_{\mathrm{PBL}}{)}_{\text{G-M}}$ ratio in Eq. (3) is provided by GEM-MACH, and has hourly temporal resolution. This modelled $(C/V_{\mathrm{PBL}}{)}_{\text{G-M}}$ ratio is referred to here as a conversion ratio $R_{\mathrm{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 ${\mathrm{NO}}_{\mathrm{2}}$ 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${}^{-\mathrm{1}}$ (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 ${\mathrm{NO}}_{\mathrm{2}}$ in the PBL) corresponds to 5–8 pptv of surface ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ concentration is controlled not only by PBL height but also by both boundary layer dynamics and photochemistry. The surface ${\mathrm{NO}}_{\mathrm{2}}$ derived using the hourly modelled $R_{\mathrm{CV}}$ ratio is referred to here as $C_{\text{pan-model}}$, while the surface ${\mathrm{NO}}_{\mathrm{2}}$ derived using the monthly mean PSC-LUT is referred to here as $C_{\text{pan-LUT}}$. In general, $C_{\text{pan-model}}$ is a data product that depends on daily model outputs, but $C_{\text{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_{\text{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_{\text{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 ${\mathrm{NO}}_{\mathrm{2}}$ data products are discussed in the next section.

Figure 5

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

[Figure omitted. See PDF]

Figure 6

Modelled and Pandora zenith-sky surface ${\mathrm{NO}}_{\mathrm{2}}$ vs. in situ ${\mathrm{NO}}_{\mathrm{2}}$ (2016–2017). Panel (a) shows the GEM-MACH modelled surface ${\mathrm{NO}}_{\mathrm{2}}$ data vs. in situ ${\mathrm{NO}}_{\mathrm{2}}$; panels (b) and (c) show the Pandora ZS surface ${\mathrm{NO}}_{\mathrm{2}}$ data vs. in situ ${\mathrm{NO}}_{\mathrm{2}}$. The Pandora ZS surface ${\mathrm{NO}}_{\mathrm{2}}$ 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.

[Figure omitted. See PDF]

4.2 Comparison with measurements and model

Figure 6 shows the evaluation of modelled and Pandora zenith-sky surface ${\mathrm{NO}}_{\mathrm{2}}$ concentrations, both using in situ ${\mathrm{NO}}_{\mathrm{2}}$ 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=\mathrm{0.78}$) and moderate positive bias (37 %, Fig. 6a). The Pandora zenith-sky surface ${\mathrm{NO}}_{\mathrm{2}}$ data, $C_{\text{pan-model}}$, shows almost the same correlation ($R=\mathrm{0.77}$), with only $-\mathrm{7}$ % bias (Fig. 6b). The better performance of $C_{\text{pan-model}}$ is expected since the conversion method for Pandora zenith-sky measurements relies on the GEM-MACH modelled ${\mathrm{NO}}_{\mathrm{2}}$ profile (see Eq. 3); in other words, the Pandora zenith-sky surface ${\mathrm{NO}}_{\mathrm{2}}$ has at least one more piece of information (i.e., ${\mathrm{NO}}_{\mathrm{2}}$ total column) than GEM-MACH surface ${\mathrm{NO}}_{\mathrm{2}}$ concentrations. The $C_{\text{pan-LUT}}$ shows a similar correlation coefficient ($R=\mathrm{0.73}$) and has improved bias ($-\mathrm{3}$ %, Fig. 6c). This result (slightly lower correlation) is also reasonable and acceptable since $C_{\text{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 ${\mathrm{NO}}_{\mathrm{2}}$ concentrations, $C_{\text{pan-model}}$ and $C_{\text{pan-LUT}}$ (Fig. 6e and f), also have better frequency distributions than the GEM-MACH (Fig. 6d). Figure 6d shows that the ${\mathrm{NO}}_{\mathrm{2}}$ surface concentrations peaks (ambient background concentrations) from model and in situ data are misaligned. This indicates that the GEM-MACH ${\mathrm{NO}}_{\mathrm{2}}$ background surface concentrations have a 1 ppbv low bias at this site. In contrast, the zenith-sky surface ${\mathrm{NO}}_{\mathrm{2}}$ at peak frequency matches the in situ data (Fig. 6e and f), indicating that the low bias of the background surface ${\mathrm{NO}}_{\mathrm{2}}$ value has been corrected with this additional information from Pandora zenith-sky total column measurements. In addition, in high ${\mathrm{NO}}_{\mathrm{2}}$ concentration conditions (> 20 ppbv), the zenith-sky surface ${\mathrm{NO}}_{\mathrm{2}}$ also shows better agreement with the in situ ${\mathrm{NO}}_{\mathrm{2}}$ than do the modelled data. The mean of the top 10 % of the in situ data is $\mathrm{26}\pm \mathrm{1}$ ppbv (uncertainty of the mean), whereas the corresponding values for GEM-MACH, $C_{\text{pan-model}}$, and $C_{\text{pan-LUT}}$ are $\mathrm{39}\pm \mathrm{1}$, $\mathrm{26}\pm \mathrm{1}$, and $\mathrm{27}\pm \mathrm{1}$ ppbv, respectively.

Figure 7

Modelled and Pandora direct-Sun surface ${\mathrm{NO}}_{\mathrm{2}}$ vs. in situ ${\mathrm{NO}}_{\mathrm{2}}$ (2016–2017). Panel (a) shows the GEM-MACH modelled surface ${\mathrm{NO}}_{\mathrm{2}}$ data vs. in situ ${\mathrm{NO}}_{\mathrm{2}}$; panels (b) and (c) show the Pandora DS surface ${\mathrm{NO}}_{\mathrm{2}}$ data vs. in situ ${\mathrm{NO}}_{\mathrm{2}}$. The Pandora DS surface ${\mathrm{NO}}_{\mathrm{2}}$ 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.

[Figure omitted. See PDF]

The total column-to-surface concentration conversion algorithm has also been applied to the Pandora direct-Sun total column ${\mathrm{NO}}_{\mathrm{2}}$ (see Fig. 7). Figure 7b shows that the direct-Sun surface ${\mathrm{NO}}_{\mathrm{2}}$ data have a similar agreement with the in situ data ($-\mathrm{8}$ % bias and $R=\mathrm{0.80}$) as the zenith-sky surface ${\mathrm{NO}}_{\mathrm{2}}$. In high ${\mathrm{NO}}_{\mathrm{2}}$ 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 $\mathrm{27}\pm \mathrm{1}$ ppbv, whereas the corresponding values for GEM-MACH, $C_{\text{pan-model}}$, and $C_{\text{pan-LUT}}$ are $\mathrm{40}\pm \mathrm{1}$, $\mathrm{27}\pm \mathrm{1}$, and $\mathrm{27}\pm \mathrm{1}$ ppbv, respectively.

Thus, in general, both Pandora zenith-sky and direct-Sun surface ${\mathrm{NO}}_{\mathrm{2}}$ datasets can be used reliably to obtain surface concentrations. The good consistency between $C_{\text{pan-model}}$ and $C_{\text{pan-LUT}}$ implies that two versions of Pandora surface ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ in heavy cloud conditions. Due to enhanced scattering, heavy clouds could lead to a significant overestimation of surface ${\mathrm{NO}}_{\mathrm{2}}$ derived from zenith-sky measurements. A cloud filtering method based on retrieved ${\mathrm{O}}_{\mathrm{4}}$ dSCDs is used to identify these conditions. High retrieved ${\mathrm{O}}_{\mathrm{4}}$ values correspond to long optical path lengths, and therefore it is expected that corresponding ${\mathrm{NO}}_{\mathrm{2}}$ values are overestimated as discussed in Appendix C.

Figure 8

Example of surface ${\mathrm{NO}}_{\mathrm{2}}$ concentration time series in all conditions (April 2017). The in situ, Pandora DS, and Pandora ZS surface ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{O}}_{\mathrm{4}}$ (heavy cloud indicator) are also labelled by green squares. Dates are in mm/dd format.

[Figure omitted. See PDF]

The effectiveness of the zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{O}}_{\mathrm{4}}$ (Appendix C), Pandora zenith-sky-derived surface ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ data can both be produced. However, when Sun strength is moderate (30–60), only zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ data are reliable. When Sun strength is low (< 30), zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ profile (this factor only applies to Pandora surface ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ profiles, surface ${\mathrm{NO}}_{\mathrm{2}}$ concentrations were obtained that agree well with in situ measurements collected at the same location. The Pandora surface ${\mathrm{NO}}_{\mathrm{2}}$ data were also analyzed in more detail with a focus on temporal variations. Figure 9 shows the averaged surface ${\mathrm{NO}}_{\mathrm{2}}$ diurnal variations of four different datasets. The in situ instrument produces continuous measurements 24 h d${}^{-\mathrm{1}}$, whereas Pandora only has measurements when sunlight is available. The diurnal variation of surface ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 $\mathrm{1}\mathit{\sigma }$ envelope.

[Figure omitted. See PDF]

Figure 9 shows that all four datasets/curves captured the enhanced morning surface ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ is $\mathrm{14.9}\pm \mathrm{9.3}$ ppbv, while GEM-MACH, Pandora DS, and Pandora ZS ${\mathrm{NO}}_{\mathrm{2}}$ are $\mathrm{23.5}\pm \mathrm{15.0}$, $\mathrm{15.6}\pm \mathrm{10.5}$, and $\mathrm{15.2}\pm \mathrm{6.8}$ ppbv, respectively. At 17:00 LST, in situ ${\mathrm{NO}}_{\mathrm{2}}$ is $\mathrm{7.3}\pm \mathrm{5.8}$ ppbv, while GEM-MACH, Pandora DS, and Pandora ZS ${\mathrm{NO}}_{\mathrm{2}}$ are $\mathrm{5.6}\pm \mathrm{5.0}$, $\mathrm{3.6}\pm \mathrm{2.6}$, and $\mathrm{5.2}\pm \mathrm{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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 $\mathrm{1}\mathit{\sigma }$ envelope.

[Figure omitted. See PDF]

The diurnal variation is also examined by grouping the data by seasons. Figure 10 shows that the surface ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ data show good agreement with in situ measurements in all seasons; the hourly mean values of Pandora surface ${\mathrm{NO}}_{\mathrm{2}}$ are all well within the $\mathrm{1}\mathit{\sigma }$ envelope of the in situ measurements.

5.2 Planetary boundary layer effect

The larger morning offset in modelled surface ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$. Figure 11a shows that, in general, the difference between modelled and in situ ${\mathrm{NO}}_{\mathrm{2}}$ decreases with an increase of PBL height. When the modelled PBL height is less than 100 m, the mean difference is $\mathrm{18}\pm \mathrm{12}$ ppbv ($\mathrm{1}\mathit{\sigma }$), while when the modelled PBL height is 1 km, the mean difference is only $\mathrm{2.9}\pm \mathrm{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 ${\mathrm{NO}}_{\mathrm{2}}$ concentrations; the $x$ axes for the right column are the difference between Pandora zenith-sky ($C_{\text{pan-model}}$) and in situ surface ${\mathrm{NO}}_{\mathrm{2}}$ 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).

[Figure omitted. See PDF]

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 ${\mathrm{NO}}_{\mathrm{2}}$ concentrations derived from Pandora zenith-sky data ($C_{\text{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 $\mathrm{0.9}\pm \mathrm{8.9}$ ppbv. When the modelled PBL height is 1 km, the mean difference is slightly improved to $\mathrm{0.1}\pm \mathrm{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 ${\mathrm{NO}}_{\mathrm{2}}$ has larger bias compared to in situ data in the morning. Figure 11d, f, and h show that Pandora zenith-sky surface ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$, 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 ${\mathrm{NO}}_{\mathrm{2}}$ from direct-Sun spectral measurements in the UV–visible spectrum. In this work, a new zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ data shows only $-\mathrm{4}$ % bias compared to the standard Pandora direct-Sun data product. In addition, OMI ${\mathrm{NO}}_{\mathrm{2}}$ SPv3.0 data demonstrate similar biases ($-\mathrm{30}$ % and $-\mathrm{29}$ %, respectively) when compared to direct-Sun and zenith-sky Pandora total column ${\mathrm{NO}}_{\mathrm{2}}$ data.

Surface ${\mathrm{NO}}_{\mathrm{2}}$ concentrations were calculated from Pandora direct-Sun and zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$ using column-to-surface ratios derived from GEM-MACH regional chemical transport model. The bias between Pandora-based direct-Sun and zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ surface concentration estimates and in situ measurements is only $-\mathrm{8}$ % and $-\mathrm{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 ${\mathrm{NO}}_{\mathrm{2}}$ concentrations also show good diurnal and seasonal variation when compared to the in situ data. High surface ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ in the morning by 8.6 ppbv (at 07:00 LST). It appears that the bias in modelled surface ${\mathrm{NO}}_{\mathrm{2}}$ is related at least in part to an incorrectly diagnosed PBL height. In contrast, the difference between Pandora-based and in situ ${\mathrm{NO}}_{\mathrm{2}}$ does not show any significant dependence on the PBL height. Thus, to enable a fast and practical Pandora surface ${\mathrm{NO}}_{\mathrm{2}}$ data production, the use of a pre-calculated conversion ratio PSC-LUT is recommended.

The new retrieval algorithm for Pandora zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ measurements can provide high-quality ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ concentration data can be obtained for air quality monitoring purposes. The variation of conversion ratios in the PSC-LUT demonstrates that the surface ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{O}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}$ ratio algorithm (Cede, 2019), it may be possible for Pandora measurements to be used to derive ${\mathrm{NO}}_{\mathrm{2}}$ tropospheric profiles and columns.

Data availability

Pandora data are available from the Pandonia network (http://pandonia.net/media/documents/BlickSoftwareSuite_Manual_v11.pdf, Cede, 2019). In situ surface ${\mathrm{NO}}_{\mathrm{2}}$ data are available from the National Air Pollution Surveillance (NAPS) program (http://maps-cartes.ec.gc.ca/rnspa-naps/data.aspx, last access: 15 August 2019). OMI ${\mathrm{NO}}_{\mathrm{2}}$ SPv3.0 data are available from 10.5067/Aura/OMI/DATA2017 (Krotkov et al., 2019). Any additional data may be obtained from Xiaoyi Zhao ([email protected]).

Appendix A Empirical zenith-sky AMF

Before calculating the empirical zenith-sky AMF, the VCD${}_{\mathrm{DS}}$ and dSCD${}_{\mathrm{ZS}}$ have both been strictly filtered to ensure any measurements used in this calculation have the highest quality. For VCD${}_{\mathrm{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${}_{\mathrm{ZS}}$, data are filtered using similar criteria as for VCD${}_{\mathrm{DS}}$, with adjustments for zenith-sky observations.

The VCD${}_{\mathrm{DS}}$ and dSCD${}_{\mathrm{ZS}}$ data are merged and divided into several SZA bins. Each bin covers 5${}^{\circ }$. A multi-non-linear regression is performed by using the following equation:

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}$$ where VCD${}_n$ is not a single direct-Sun VCD data point but is an $m\times \mathrm{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${}^{\circ }$ SZA bin, and each element of the $m\times \mathrm{1}$ matrix is a single VCD in that SZA bin. Similarly, dSCD${}_n$ is also an $m\times \mathrm{1}$ matrix, with each element representing a single coincident zenith-sky dSCD in SZA bin number $n$. $I_n$ is an $m\times \mathrm{1}$ indicator function, where the elements of $I_n$ are set to 1. The RCD and $b_{\mathrm{1}}$ 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_{\mathrm{1}}$ 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${}^{\circ }$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ AMF function, which has the expression A2$$\mathrm{AMF}=a_{\mathrm{1}}+\left(\mathrm{1.02}-a_{\mathrm{1}}\right)/\cos \left(\mathrm{SZA}\right).$$ 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${}^{\circ }$). 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 ${\mathrm{NO}}_{\mathrm{2}}$ has to be accounted (i.e., the geometry-form AMFs are not enough).

Appendix B

Stratospheric ${\mathrm{NO}}_{\mathrm{2}}$ column

Several stratospheric ${\mathrm{NO}}_{\mathrm{2}}$ column values were tested and used in the surface ${\mathrm{NO}}_{\mathrm{2}}$ concentration algorithm (Eq. 3). Figure B1a shows the OMI monthly mean (referred to as OMI) and Pratmo box-model stratospheric column ${\mathrm{NO}}_{\mathrm{2}}$ (Adams et al., 2016; McLinden et al., 2000) (referred to as box). Since the satellite only samples Toronto once per day, the OMI stratospheric ${\mathrm{NO}}_{\mathrm{2}}$ lacks diurnal variation. To account for the diurnal variation, diurnal ratios of ${\mathrm{NO}}_{\mathrm{2}}$ VCD have been calculated and applied to OMI monthly mean data. The stratospheric ${\mathrm{NO}}_{\mathrm{2}}$ columns are calculated using 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),$$ where $V_{\mathrm{OMI}}(t_{\mathrm{0}}$) is the OMI-measured stratospheric column, $t_{\mathrm{0}}$ is OMI overpass time, $V_{\mathrm{box}}\left(t_{\mathrm{0}}\right)$ is the modelled stratospheric column at OMI overpass time, $V_{\mathrm{box}}\left(t\right)$ is the modelled stratospheric column at time $t$, and $V_{\mathrm{OMI}}\left(t\right)$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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.

[Figure omitted. See PDF]

To justify why this diurnal variation has to be included, Fig. B1c shows the total column ${\mathrm{NO}}_{\mathrm{2}}$ time series. The diurnal stratospheric ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ data (e.g., in the morning, this will lead to the overestimation of the stratospheric ${\mathrm{NO}}_{\mathrm{2}}$ and thus the underestimation of surface ${\mathrm{NO}}_{\mathrm{2}}$). Please note that the strength of this bias is related to (1) the ${\mathrm{NO}}_{\mathrm{2}}$ profile (weights between stratospheric and tropospheric ${\mathrm{NO}}_{\mathrm{2}}$), 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$. Depending on the properties of the clouds and the ${\mathrm{NO}}_{\mathrm{2}}$ profile, the clouds could have non-negligible impacts on zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ observations.

A typical method of removing zenith-sky measurements affected by heavy clouds is to eliminate measurements with large enhancements of ${\mathrm{O}}_{\mathrm{4}}$ and/or ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ (Van Roozendael and Hendrick, 2012). In the Pandora zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ retrieval, we use the ${\mathrm{O}}_{\mathrm{4}}$ dSCDs. Since the measured ${\mathrm{O}}_{\mathrm{4}}$ dSCDs has SZA dependency, all measured ${\mathrm{O}}_{\mathrm{4}}$ 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 ${\mathrm{O}}_{\mathrm{4}}$ dSCDs.

Figure C1

Illustration of cloud effect and heavy cloud data filtration: panel (a) shows measured ${\mathrm{O}}_{\mathrm{4}}$ differential slant column densities vs. solar zenith angle; the grey dots represent the top 0–10th percentile range of ${\mathrm{O}}_{\mathrm{4}}$. Panel (b) shows the scatter plot of zenith-sky vs. in situ surface ${\mathrm{NO}}_{\mathrm{2}}$ data that has ${\mathrm{O}}_{\mathrm{4}}$ 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 ${\mathrm{O}}_{\mathrm{4}}$; panel (d) is similar to (b) but uses the data that has ${\mathrm{O}}_{\mathrm{4}}$ 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 omitted. See PDF]

Figure C2

Correlation coefficient and bias (slope) between zenith-sky and in situ surface ${\mathrm{NO}}_{\mathrm{2}}$ data in different ${\mathrm{O}}_{\mathrm{4}}$ dSCD percentile bins. Panel (a) shows the correlation coefficients; panel (b) shows the slopes of linear fit with intercept set to 0.

[Figure omitted. See PDF]

Figure C3

Example of surface ${\mathrm{NO}}_{\mathrm{2}}$ concentration time series in all conditions. The in situ, Pandora DS, and Pandora ZS surface ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{O}}_{\mathrm{4}}$ (heavy cloud indicator) are also labelled by green squares. Dates are in mm/dd format.

[Figure omitted. See PDF]

Figure C1 shows examples of selected Pandora zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ data and their corresponding ${\mathrm{O}}_{\mathrm{4}}$ dSCD values. For example, Fig. C1a shows the ${\mathrm{O}}_{\mathrm{4}}$ dSCDs vs. SZA, and the top 10 percentiles of the data with enhanced ${\mathrm{O}}_{\mathrm{4}}$ are marked in grey. The corresponding Pandora zenith-sky data are plotted against in situ data in Fig. C1b, which shows low correlation ($R=\mathrm{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 ${\mathrm{NO}}_{\mathrm{2}}$ retrieval. Figures C1c and d are similar to Fig. C1a and b but for selected Pandora zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ data that have ${\mathrm{O}}_{\mathrm{4}}$ values within the 40th to 50th percentile range. Figure C1d shows that when ${\mathrm{O}}_{\mathrm{4}}$ is not enhanced, the derived zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$ has good agreement with in situ data ($R=\mathrm{0.8}$ and bias $=-\mathrm{5}$ %). To summarize how the retrieved ${\mathrm{O}}_{\mathrm{4}}$ dSCDs can indicate the quality of the Pandora zenith-sky ${\mathrm{NO}}_{\mathrm{2}}$, 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 ${\mathrm{NO}}_{\mathrm{2}}$ data that have ${\mathrm{O}}_{\mathrm{4}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 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}}},$$ where ${\mathit{\sigma }}_{\mathrm{dSCD}}$ is the statistical uncertainty on the DOAS fit (output of QDOAS) and ${\mathit{\sigma }}_{\mathrm{RCD}}$ and ${\mathit{\sigma }}_{\mathrm{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 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}$$ where ${\mathit{\sigma }}_{V_{\mathrm{pan}}}$ is the uncertainty of Pandora zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$, (here we use the derived ${\mathit{\sigma }}_{\mathrm{VCD}}$ in Eq. D1), ${\mathit{\sigma }}_{V_{\mathrm{strat}}}$ is the uncertainty of the stratospheric ${\mathrm{NO}}_{\mathrm{2}}$ column (estimated using the $\mathrm{1}\mathit{\sigma }$ standard deviation of the $V_{\mathrm{strat}}$), ${\mathit{\sigma }}_{V_{\mathrm{ftrop}}}$ is the uncertainty of the free troposphere ${\mathrm{NO}}_{\mathrm{2}}$ column (estimated using the $\mathrm{1}\mathit{\sigma }$ standard deviation of the $V_{\mathrm{ftrop}}$). $R_{\mathrm{CV}}$ is the GEM-MACH calculated surface VMR to PBL column ratio, and ${\mathit{\sigma }}_{\mathrm{R}}$ is the uncertainty of that ratio (estimated using the $\mathrm{1}\mathit{\sigma }$ standard deviation of the $R_{\mathrm{CV}}$). The means of $R_{\mathrm{CV}}$, $V_{\mathrm{Pan}}$, $V_{\mathrm{strat}}$, and $V_{\mathrm{ftrop}}$ are used in the uncertainty estimation.

Table D1

Estimated uncertainties for Pandora zenith-sky total column and surface ${\mathrm{NO}}_{\mathrm{2}}$.

Estimation	${\mathit{\sigma }}_{{\mathrm{VCD}}_{\mathrm{ZS}}}$	${\mathit{\sigma }}_{C_{\text{Pan-ZS}}}$	${\mathit{\sigma }}_{C_{\text{Pan-DS}}}$
method	(DU)	(ppbv)	(ppbv)
UP	0.12	6.5	5.4
SDD	0.09	5.1	5.0
Unbiased SDD	0.09	4.8	4.8

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 ${\mathrm{NO}}_{\mathrm{2}}$ data can represent the true value, we can estimate the uncertainty of Pandora zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$ by calculating the $\mathrm{1}\mathit{\sigma }$ standard deviation of their difference (referred to here as the SDD method): D3$${\mathit{\sigma }}_{{\mathrm{VCD}}_{\mathrm{ZS}}}=\mathit{\sigma }\left({\mathrm{VCD}}_{\mathrm{DS}}-{\mathrm{VCD}}_{\mathrm{ZS}}\right).$$ Similarly, if we assume that the in situ surface ${\mathrm{NO}}_{\mathrm{2}}$ VMR can represent the true value, the uncertainty of Pandora zenith-sky-based surface ${\mathrm{NO}}_{\mathrm{2}}$ VMR can be given by D4$${\mathit{\sigma }}_{C_{\mathrm{Pan}}}=\mathit{\sigma }\left(C_{\text{in situ}}-C_{\mathrm{pan}}\right).$$ 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_{\mathrm{1}}$ and $k_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ 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 $\mathrm{1}\mathit{\sigma }$ level). When using the UP method, for the worst-case scenario, the Pandora zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$ 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 ${\mathrm{NO}}_{\mathrm{2}}$ VMR data have uncertainties from 4.8 to 6.5 ppbv. In Eq. (D2), the contributions of the $V_{\mathrm{Pan}}$, $V_{\mathrm{Strat}}$, $V_{\mathrm{ftrop}}$, and $R_{\mathrm{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 ${\mathrm{NO}}_{\mathrm{2}}$ VMR is dominated by the uncertainties of Pandora zenith-sky total column ${\mathrm{NO}}_{\mathrm{2}}$ and the modelled column-to-surface conversion ratio ($R_{\mathrm{CV}}$). However, note that this uncertainty budget depends on the ${\mathrm{NO}}_{\mathrm{2}}$ vertical distributions and hence may vary from site to site; e.g., in Toronto, tropospheric column ${\mathrm{NO}}_{\mathrm{2}}$ is typically 2–4 times higher than stratospheric column ${\mathrm{NO}}_{\mathrm{2}}$, and thus the contribution to uncertainty from $V_{\mathrm{Pan}}$ is much larger than the corresponding contributions from $V_{\mathrm{Strat}}$ and $V_{\mathrm{ftrop}}$. In addition, the uncertainty of Pandora direct-Sun surface ${\mathrm{NO}}_{\mathrm{2}}$ VMR is also estimated and provided in Table D1. It shows slightly better results than for zenith-sky-based surface ${\mathrm{NO}}_{\mathrm{2}}$ VMR.

Author contributions

XZ analyzed the data and prepared the manuscript, with significant conceptual input from DG, VF, and CM, and critical feedback from all co-authors. JD, AO, VF, XZ, and SCL operated and managed the Canadian Pandora network. AL, MDM, and DG performed and analyzed the GEM-MACH simulations. AC, MT, and MM operated the Pandonia network and provided critical technical support to the Canadian Pandora measurements and subsequent data analysis.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

Xiaoyi Zhao was supported by the NSERC Visiting Fellowships in Canadian Government Laboratories program. We thank Ihab Abboud and Reno Sit for their technical support of Pandora measurements. We thank NAPS for providing surface ${\mathrm{NO}}_{\mathrm{2}}$ data. We acknowledge the NASA Earth Science Division for providing OMI ${\mathrm{NO}}_{\mathrm{2}}$ SPv3.0 data. We also thank Thomas Danckaert, Caroline Fayt, Michel Van Roozendael, and others from IASB-BIRA for providing the QDOAS software, the NDACC UV–visible working group for providing NDACC UV–visible ${\mathrm{NO}}_{\mathrm{2}}$ AMF LUT, and Yushan Su from the Ontario Ministry of the Environment, Conservation and Parks for providing NAPS Toronto north station in situ ${\mathrm{NO}}_{\mathrm{2}}$ information. We thank two anonymous referees for their helpful and insightful comments, which improved the overall quality of this work.

Review statement

This paper was edited by Robert McLaren and reviewed by Michel Van Roozendael and one anonymous referee.