Journal topic
Atmos. Meas. Tech., 12, 2545–2565, 2019
https://doi.org/10.5194/amt-12-2545-2019
Atmos. Meas. Tech., 12, 2545–2565, 2019
https://doi.org/10.5194/amt-12-2545-2019

Research article 02 May 2019

Research article | 02 May 2019

# NO2 vertical profiles and column densities from MAX-DOAS measurements in Mexico City

NO2 vertical profiles and column densities from MAX-DOAS measurements in Mexico City
Martina Michaela Friedrich1,2, Claudia Rivera1,3, Wolfgang Stremme1, Zuleica Ojeda1, Josué Arellano1, Alejandro Bezanilla1, José Agustín García-Reynoso1, and Michel Grutter1 Martina Michaela Friedrich et al.
• 1Centro de Ciencias de la Atmósfera, Universidad Nacional Autónoma de México, Mexico City, Mexico
• 2Belgian Institute for Space Aeronomie (BIRA-IASB), Brussels, Belgium
• 3Facultad de Química, Universidad Nacional Autónoma de México, Mexico City, Mexico

Correspondence: Martina Michaela Friedrich (martina@atmosfera.unam.mx)

Abstract

We present a new numerical code, Mexican MAX-DOAS Fit (MMF), developed to retrieve profiles of different trace gases from the network of MAX-DOAS instruments operated in Mexico City. MMF uses differential slant column densities (dSCDs) retrieved with the QDOAS software. The retrieval is comprised of two steps, an aerosol retrieval and a trace gas retrieval that uses the retrieved aerosol profile in the forward model for the trace gas. For forward model simulations, VLIDORT is used (e.g., Spurr et al.2001; Spurr2006, 2013). Both steps use constrained least-square fitting, but the aerosol retrieval uses Tikhonov regularization and the trace gas retrieval optimal estimation. Aerosol optical depth and scattering properties from the AERONET database, averaged ceilometer data, WRF-Chem model data, and temperature and pressure sounding data are used for different steps in the retrieval chain.

The MMF code was applied to retrieve NO2 profiles with 2 degrees of freedom (DOF = 2) from spectra of the MAX-DOAS instrument located at the Universidad Nacional Autónoma de México (UNAM) campus. We describe the full error analysis of the retrievals and include a sensitivity exercise to quantify the contribution of the uncertainties in the aerosol extinction profiles to the total error. A data set comprised of measurements from January 2015 to July 2016 was processed and the results compared to independent surface measurements. We concentrate on the analysis of four single days and additionally present diurnal and annual variabilities from averaging the 1.5 years of data. The total error, depending on the exact counting, is 14 %–20 % and this work provides new and relevant information about NO2 in the boundary layer of Mexico City.

1 Introduction

Air pollution is a serious environmental problem due to its negative impacts on human health and ecosystems. Fast-growing urban and industrial centers are continuously affected by bad air quality, and in order to assess their current efforts to mitigate emissions and plan for more efficient strategies to lower the concentration levels of harmful contaminants, it is indispensable to have the proper tools to measure them not only at ground level but also throughout the boundary layer . The multi-axis differential optical absorption spectroscopy (MAX-DOAS) technique (e.g., Hönninger et al.2004; Platt and Stutz2008) has rapidly developed in recent years and has proven extremely valuable in tropospheric chemistry and air pollution studies since it provides vertical distribution of trace gases with high temporal resolution.

Figure 1Flowchart of the complete trace gas retrieval algorithm. Orange boxes belong to the aerosol retrieval and green boxes to trace gas retrieval. The light blue box encompasses forward model input calculation. The dark blue boxes are in- and outputs of QDOAS. The yellow boxes represent the forward modeling steps. The red boxes are the inversion steps, using Tikhonov regularization for aerosol retrieval and optical estimation (OE) for trace gas retrieval.

This remote-sensing technique is based on the spectroscopically resolved measurement of scattered sunlight at different elevation angles, allowing for the retrieval of total column amounts of aerosols and trace gases with profiling capability. Powerful applications of this technique have been demonstrated to provide useful information about the vertical distribution of aerosols and photochemically relevant species such as nitrogen dioxide (NO2), formaldehyde (HCHO), glyoxal (CHOCHO) and nitrous acid (HONO) among other gases .

Photochemical reactions involving NO2 play an important role in the formation of O3 . The Mexico City metropolitan area (MCMA) has been particularly affected since the 1990s by high-O3 episodes threatening the population and forcing the authorities to impose strict restrictions on the usage of motor vehicles . Measurements have been performed in the region using fixed and mobile DOAS zenith-scattered sunlight . In 2014, a MAX-DOAS network, initially consisting of four instruments, was established in the Mexico City metropolitan area (MCMA) and has been operating since .

In this contribution, we describe the MMF (Mexican MAX-DOAS Fit) code that has been implemented to retrieve vertical distribution of aerosols and trace gases with emphasis on the errors and diagnostics of the results. An overview of the MAX-DOAS instruments is provided in Sect. 2 and the complete retrieval strategy from the measured spectra to vertical trace gas profiles is summarized in Fig. 1.

Radiative transfer simulations (constituting the forward model in MMF, yellow boxes in Fig. 1) are performed with VLIDORT (Spurr2013) to derive simulated differential slant column densities (dSCDs) at the middle of the corresponding wavelength interval used to derive dSCDs from measured spectra. The dSCD retrieval (blue boxes in Fig. 1) is performed with QDOAS and is described briefly in Sect. 3. The orange parts in the figure refer to the aerosol retrieval while the green parts belong to the trace gas retrieval. Details on the forward model choice, the forward model input calculation (light blue box in Fig. 1) and processing of output quantities in the inversion algorithm are described in Sect. 4.

An error analysis has been included in Sect. 5, in particular investigating the effect of the aerosol retrieval on the NO2 results. Some examples of the NO2 variability are provided from one of the network's stations and compared to surface concentrations in Sect. 6, and a summary of the work and an outlook on planned improvements to the retrieval code are presented in Sect. 7.

2 Instruments

An instrument based on the MAX-DOAS technique was designed and developed by the Center for Atmospheric Sciences at the Universidad Nacional Autónoma de México (UNAM). It consists of two main parts: the scanner unit, which collects the scattered light, and the acquisition–control unit, containing a spectrometer and computer that record and store the measurements. The two components are connected by an optical fiber and a data connection cable. Both are described briefly in the following sections. For a more complete description and a map of the area with the station locations and their measurement orientation, we refer the reader to .

## 2.1 The scanner

The scanner is composed of a plastic enclosure (NEMA-rated type 3) resistant to sunlight and hermetically sealed to protect the inner parts from water and bugs. This is important in order to assure long measurement periods under harsh conditions . The scattered light is collected by a plano-convex quartz lens (Edmund Optics, =25.4 mm, f=100 mm), focusing it into the entrance of an optical fiber (Fibertech Optica, quartz, 6 m long, =0.6 mm), which transfers the light into the acquisition–control unit. These optical parts are mounted inside a telescope housing constructed of Nylamid material. A shutter system is also installed within the telescope. It consists of a stepper motor (Mercury, 7.5 by step) that rotates a metal circular plate to prevent the passage of light into the optical fiber and an optical switch which is used to indicate the position of the plate. The shutter is used to make measurements of dark spectra in between scans.

A stepper motor (Oriental Motors, PK266-02A) inside the enclosure allows the movement of the scanner unit in a range of 180 with steps of 0.1. A mechanical switch is used as a reference to indicate the starting position of each scan. The motor is controlled by an electronic board composed of an 8-bit microcontroller (AVR architecture), a RS-232 port for communication with the acquisition–control unit, a temperature sensor (Maxim 18B20, accuracy ±0.5 C) and a dual-axis accelerometer (Analog Devices, accuracy ±0.1) that provides an accurate determination of the telescope's pointing elevation. The theoretical field of view (FOV) of the optical system is 0.31.

## 2.2 The acquisition–control unit

The second part of the instrument consists of a metallic housing receiving the collected light through the optical fiber and sending it to the spectrometer (Ocean Optics, USB2000+). This commercial device has a crossed asymmetric Czerny–Turner configuration, diffraction grating (1800 lines mm−1) and a slit size of 50 µm wide × 1 mm high, recording spectra in a wavelength range of 289–510 nm at a resolution of 0.69 nm (full width at half maximum). It uses a charge-coupled device (CCD) detector array (Sony ILX511B) of 2048 pixels with an integration time adjustable between 1 ms and 65 s.

Because changes in temperature can affect the wavelength  pixel ratio and also modify the optical properties of the spectrometer like the alignment and the line shape (e.g., Carlson et al.2010; Coburn et al.2011), a temperature control system composed of a Peltier cell (Multicomp) and three temperature sensors (Maxim 18B20, accuracy ±0.5C) controlled by an electronic board was implemented. The cooling side of the Peltier cell was placed on top of the spectrometer and the heating side was attached to a heat sink. The Peltier cell and the spectrometer were wrapped in a styrofoam box to keep the temperature insulated from the outside. The three temperature sensors were placed on the heat sink, the Peltier cell and the spectrometer to monitor temperature changes. The temperature control system was wrapped in an aluminum enclosure to prevent the heat spreading to other parts of the system and a fan was installed to extract the heat from the enclosure.

The electronic board for the temperature control is composed of an 8-bit microcontroller (AVR, architecture) and a RS-232 communication port. This board is responsible for obtaining the data from the sensors and adjusting the voltage in the Peltier cell to keep the temperature constant.

The acquisition–control unit has a laptop computer (Dell, Latitude 2021) with a Linux operating system contained within the enclosure. The program controlling the hardware is written in C++ using Qt libraries. A script is used to carry out the measurement sequence previously defined and to monitor the spectrometer temperature.

## 2.3 Measurement strategy

A complete scan consists of a sequence that begins with a measurement towards the zenith (90 elevation angle), followed by tree measurements between elevation angles of 0 and 10 towards the west (for UNAM station, this is an 85 azimuth angle); then 12 measurements are taken with elevation angles between 10 towards the west and 10 towards the east (crossing the zenith, but without taking a measurement), followed by three measurements between 10 and 0 elevation angle towards the east. The same sequence is then repeated but in reverse order. At the end of this cycle, a dark spectrum with the closed shutter is taken. With this setup, a complete scan takes about 7 min. The measurement sequence (90 zenith, 0, 2, 6, 13, 23, 36, 50, 65, 82 W, 82, 65, 50, 36, 23, 13, 6, 2, 0 E) is likely to change in the future to use longer integration times but at the same time reduce the number of viewing directions in the range between a 10 elevation angle towards the west and a 10 elevation angle towards the east in order to keep the total scan time roughly constant.

All data presented in this paper use this setup, meaning that we include both westerly and easterly viewing directions in the same retrieval and hence the result represents an average. This strategy leads to larger fitting errors since this likely deviates further from the assumption of horizontal homogeneity. With our retrieval chain, it is possible to consider the easterly and westerly directions separately to investigate the differences in viewing directions, which is the subject of current investigation.

The output from the MAX-DOAS instruments consists of five files per complete scan: a file containing the metadata of each spectrum (e.g., accelerometer data for each measurement, time of the acquisition, temperatures within the acquisition unit), all the spectra in non-zenith directions, the dark spectrum measured with the shutter closed, a metafile for the dark spectrum and the first zenith reference measurement. These files are stored for further processing.

3 Differential slant column density (dSCD) retrieval

The spectra were evaluated using the QDOAS (version 2.105) software . As a preprocessing step before the QDOAS analysis, the dark signal was subtracted from each of the measurement spectra. A wavelength calibration was conducted in QDOAS by applying a nonlinear least-squares fit to a solar atlas .

For NO2, the retrieval was conducted in the 405 to 465 nm wavelength range using the spectrum measured at the zenith position at the beginning of each of the measurement sequences as reference. Using a single zenith reference before the scan (both for NO2 and O4) instead of a linear interpolation of zenith measurements before and after the scan makes correct measurements during twilight problematic: since upper atmospheric contributions to the slant column density change rapidly during twilight, stratospheric NO2 could falsely be attributed to the troposphere. Hence, in the study presented here, no twilight measurements are included. For the analysis, differential cross sections of NO2 at 298 K , O3 at 221 and 241 K , and the oxygen dimer were convolved with the slit function of the spectrometer and a wavelength calibration file (created using a mercury lamp) and using the convolution tool of the QDOAS software. A Ring spectrum, generated at 273 K from a high-resolution Kurucz file using the Ring tool of the QDOAS software , was also included in the analysis.

For O4, the retrieval was conducted in the 336 to 390 nm wavelength range. Differential cross sections of O4 , O3 at 221 and 241 K , NO2 at 294 K , BrO at 298 K , HCHO at 298 K , and a Ring spectrum were included in the analysis. A third-degree polynomial was used for the retrievals.

Figure 2 shows the dSCD retrieval error statistics for NO2 as a function of elevation angle. The plot shows results of 31 448 scans from January to December 2016. Larger dSCD fitting errors are found for viewing angles closer to the horizon (smaller elevation angles), likely due to physical interferences during the measurements. As elevation angles approach the zenith, retrieval dSCD errors decrease considerably. The dSCD errors are the QDOAS calculated errors, a function of the fit residuals and the degree of freedom; see for details.

Figure 2The dSCD measurement error statistics for NO2 at the UNAM station as a function of elevation angle for data taken in the year 2016. The box encloses the 25–75th percentiles, the whiskers are 5 %–95 %, and the green bar is the mean and the red bar the median.

4 Mexican MAX-DOAS Fit

The method for the trace gas retrieval from slant column densities using the MMF code is comprised of two parts, an aerosol retrieval using the known O4 profile, and the trace gas retrieval (e.g., Platt and Stutz2008). Both parts consist of the same steps: a forward model and an inversion algorithm.

In Sect. 4.1, details about the inversion strategy are given. Our choice of forward model, VLIDORT (e.g., Spurr et al.2001; Spurr2006, 2013), and the input parameter calculation are detailed in Sect. 4.2. How forward model output quantities are processed for the inversion step is detailed in Sect. 4.3.

The retrieval time per aerosol and trace gas retrieval with the Mexico City setup is roughly half a minute for each scan, but highly dependent on the conditions.

## 4.1 Inversion theory

The inversion strategy relies on the fact that the problem is sufficiently linear so that in the iteration procedure, the new value for the quantity vector in question x (either the aerosol total extinction per layer or the trace gas optical depth per layer) can be calculated using a Gauss–Newton (GN) scheme1 according to Eq. (1) (Rodgers2000). This step corresponds to the red box and arrows in Fig. 1.

$\begin{array}{ll}\text{(1)}& {\mathbit{x}}_{i+\mathrm{1}}=& \phantom{\rule{0.25em}{0ex}}{\mathbit{x}}_{\text{a}}+{\left({\mathbf{K}}_{i}^{\text{T}}{\mathbf{S}}_{\mathrm{m}}^{-\mathrm{1}}{\mathbf{K}}_{i}+\mathbf{R}\right)}^{-\mathrm{1}}{\mathbf{K}}_{i}^{\text{T}}{\mathbf{S}}_{\mathrm{m}}& \left[\left(\mathbit{y}-\mathbit{F}\left({\mathbit{x}}_{i}\right)\right)-{\mathbf{K}}_{i}\left({\mathbit{x}}_{\text{a}}-{\mathbit{x}}_{i}\right)\right]\end{array}$

Here, superscript “T” denotes transposed; superscript “−1” denotes the inverse. The index i is the iteration index; the subscript “a” indicates a priori values. Sm is the measurement error covariance matrix, and y denotes the vector of measured differential slant column densities. F(xi) are the simulated differential slant column densities, calculated using the forward model with input profile xi. Both y and F(xi) are vectors of dimension, which is the number of telescope viewing angles. ${\mathbf{K}}_{i}=\partial \mathbit{F}\left(\mathbit{x}{\right)}^{l}/\partial {\mathbit{x}}_{n}$ is the Jacobian matrix at the ith iteration describing the change in simulated dSCD for viewing angle l when the profile x in layer n is varied.

In the case of optimal estimation (OE), the regularization matrix R is equal to the inverse of the a priori covariance matrix, $\mathbf{R}={\mathbf{S}}_{\text{a}}^{-\mathrm{1}}$. OE regularization is used for trace gas retrieval. Other regularization matrices are possible; see, e.g., .

For the aerosol retrieval used in this study, we use the L1 operator (R=L1TαL1) where the scaling parameter α is set to a constant value of 20 and is supplied via an input script to limit the degrees of freedom (DOFs) to just slightly above 1. Different scalings for the upper layers and lower layers could be supplied, as well as a complete regularization matrix R.

## 4.2 Forward model

Several radiative transfer codes have been developed to serve as forward models for this kind of retrieval algorithms. For example, there are Monte Carlo (MC) codes, such as AMFTRAN (Marquard1998), TRACY or PROMSAR . All these MC radiative transfer codes start with ejecting photons from the instrument and following the photon path backwards (see also, e.g., Perliski and Solomon1993; Marquard et al.2000); hence they are sometimes somewhat confusingly referred to as backward models. The advantage of MC codes is their high accuracy; the disadvantage is the rather long calculation time. SCIATRAN (e.g., Buchwitz et al.1998; Rozanov et al.2014, 2017) combines multiple-scattering radiative transfer modeling using the Picard iterative approximation and a chemistry and an inversion code. It has been heavily used in NO2 (e.g., Vidot et al.2010) and ozone (e.g., Rahpoe et al.2013b, a) retrieval from satellite observations.

Another class of radiative transfer codes, such as VLIDORT and DISORT , uses the discrete ordinate algorithm. This finite difference method is based on finding solutions for an atmosphere consisting of a number of homogeneous layers using Gaussian quadrature approximations to the integrodifferential radiative transfer equation expressed using Legendre polynomials for the phase function and Fourier series for the intensity.

The big advantage of the second class of codes is their high speed. Another advantage of VLIDORT is that it calculates not only the intensity but also analytic Jacobians. The number of calls to the program is not equal to twice the number of layers (as it would need to be for finite difference approximation), but a single call is sufficient. MMF uses the intensity part of VLIDORT, version 2.7, released in August 2014. (V)LIDORT is configured to return intensity Jacobians with respect to gas absorption layer properties or aerosol total extinction layer properties.

Figure 3Calculating the gas and air (a) and aerosol (b) contributions to the layer input parameters for VLIDORT. The green boxes are primary input to MMF, blue boxes are intermediate quantities and red boxes are final outputs that are combined for VLIDORT layer input. Here VMR denotes the a priori volume mixing ratio. Instead of a VMR, an a priori trace gas density can be supplied, too.

In MMF, the simulation grid is supplied via an input file to the code together with the corresponding a priori values of the trace gas concentration and the aerosol profile. The retrieval can be performed on the whole simulation grid or only on the lower part of the simulation grid. In this study, the simulation grid is identical to the retrieval grid and the layer distribution of the retrieval grid is not equidistant. The grid consists of 22 layers up to 25 km. The layer thickness increases from 100 m at the lowest layer to 5 km in the uppermost layer. The exact height distribution can be seen in Fig. 5. The surface albedo used in this study was set to 0.07; this value can be passed to MMF in the configuration file. However, in practice the effect for downwelling intensity and hence for the dSCD calculation was found to be negligible.

For each simulated atmospheric layer, the following input needs to be supplied to VLIDORT2:

1. total layer optical depth, τ;

2. single-scattering albedo, ω;

3. phase function expansion coefficients, β.

Since Jacobians with respect to changes of quantity χ in each layer are required, the normalized derivatives of the total primary optical quantities, τ, ω and β, with respect to this quantity need to be supplied as well. In the case of trace gas retrieval, χ=agasΔh, and in the case of aerosol retrieval, $\mathit{\chi }=\left({\mathit{\sigma }}_{\mathrm{aer}}+{a}_{\mathrm{aer}}\right)\mathrm{\Delta }h$. Here, agas and aaer are the gaseous and aerosol absorption coefficients, respectively, and σaer is the aerosol scattering coefficient. Δh is the layer thickness. Hence, for each layer the following input is additionally required:

• 4.

normalized rate of change in τ with respect to χ,

• 5.

normalized rate of change in ω with respect to χ,

• 6.

normalized rate of change in β with respect to χ.

In order to calculate this input, we first need to calculate the separate properties for the trace gases and the aerosols. Afterwards, these quantities are combined to yield the total layer quantities.

The calculation of the contribution from the trace gas and the air density (through Rayleigh scattering) for the layer inputs (one to six) is presented in Sect. 4.2.1. The aerosol contribution calculation is outlined in Sect. 4.2.2. These two sections also contain information on the source of information for the respective a priori values, as well as, for gas, the error covariance matrix.

The calculation of the contributions from gas and aerosol to yield the VLIDORT inputs (one to six) is detailed in Sect. 4.2.3.

### 4.2.1 Input for gases

In Fig. 3a, the steps to calculate the contribution of the trace gas and the contribution from Rayleigh scattering in each layer are outlined. Red boxes represent the trace gas (and air) contribution to VLIDORT input quantities, green boxes are primary inputs to the inversion code and light blue boxes are intermediate quantities that are neither direct input to MMF nor final input for the forward model.

The a priori volume mixing ratio (VMR) profile and the covariance matrix for this study are calculated from simulations covering the year 2011 using WRF-Chem v3.6. The model domain covers Mexico and surrounding seas; the domain has 200-by-100 square grid cells with approximately 28 km width and 35 vertical layers. The following WRF parameterizations were used: WRF single-moment three-class scheme (mp_physics = 3), planetary boundary layer (PBL, bl_pbl_physics = 1), and cumulus option Grell 3-D ensemble scheme (cu_physics = 5). As a land-surface option, the unified Noah land use surface model (sf_surface_physics = 2) was used. The time-varying sea surface temperature parameter was set to “on”, the grid analysis nudging parameter was set to “on” and the urban canopy model was set to “off”. The emissions used were for the RADM2 mechanism and are from the National Emissions Inventory 2008.

MMF takes pressure and temperature profiles on a separate height grid; internal interpolation to the provided retrieval grid is performed. In the processing chain implemented in Mexico, the temperature and pressure profiles from radiosonde data for the specific day are downloaded from the University of Wyoming (http://weather.uwyo.edu/upperair/sounding.html, last access: 9 April 2019) from the Mexico City International Airport, station number 76679. These temperature and pressure profiles are linearly interpolated to the simulation grid used to calculate the dry air number density ρ; see Appendix B.

The absorption cross section ξλ is taken at a wavelength in the middle of the wavelength interval used for the QDOAS retrieval (see Sect. 3): for aerosol retrieval it is 361 nm and for NO2 it is 414 nm. The same cross-section tables as for QDOAS dSCD retrievals (see Sect. 3) are used for NO2 and O4. However, since the cross section cancels out at the end of the calculation, its exact value is not too important.

There are analytical fits to the temperature dependence, using linear (e.g. Vandaele et al.2003, for NO2) or more complicated (e.g. Kirmse et al.1997) temperature dependences, for example. For a comparison study of different fitting functions, see . These could be implemented; however, the cross section is currently not temperature dependent, in agreement with the assumption made for defining dSCDs.

The trace gas absorption coefficient is calculated according to Eq. (2):

$\begin{array}{}\text{(2)}& {a}_{\mathrm{gas}}=\mathrm{VMR}\cdot {\mathit{\xi }}_{\mathit{\lambda }}\cdot \mathit{\rho }.\end{array}$

For the calculation of the depolarization ratio Δ, the main contributions are from N2, CO2, O2 and Ar. Our implementation follows . Details are given in Appendix B. For the calculation of the Rayleigh cross section QRay, we follow the implementation of (their Eq. 7.37); see also (their Eq. 19.5). The air scattering coefficient can then be calculated according to Eq. (3):

$\begin{array}{}\text{(3)}& {\mathit{\sigma }}_{\mathrm{air}}={\mathrm{Q}}_{\mathrm{Ray}}\cdot \mathit{\rho },\end{array}$

and the Rayleigh scattering expansion coefficients βl (e.g., Spurr et al.2001) can be calculated according to Eq. (4):

$\begin{array}{}\text{(4)}& {\mathit{\beta }}_{\mathrm{air}\mathrm{0}}=\mathrm{1},\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{\mathit{\beta }}_{\mathrm{air},\mathrm{1}}=\mathrm{0},\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{\mathit{\beta }}_{\mathrm{air},\mathrm{2}}=\left(\mathrm{1}-\mathrm{\Delta }\right)/\left(\mathrm{2}+\mathrm{\Delta }\right).\end{array}$

### 4.2.2 Input for aerosols

The (a priori) aerosol data for total optical depth, average single-scattering albedo ω and asymmetry parameter g (used to calculate the phase function moments) are time interpolated values from the co-located AERONET (Aerosol Robotic Network) station in Mexico City (V2, level 1.5 at http://aeronet.gsfc.nasa.gov, last access: 5 November 2018). They are also extra- and interpolated at the retrieval wavelength. The a priori shape of the profile is taken from hourly averaged ceilometer data , interpolated at the middle layer height h of each layer.

The first part of MMF, the aerosol retrieval, is limited to the aerosol density profile and the total aerosol optical depth. The average single-scattering albedo ω and asymmetry parameter g are not subject to retrieval and are assumed to be constant in all layers.

Figure 3b outlines the strategy how the aerosol contribution to the VLIDORT layer input parameters is calculated.

The hourly averages of (relative) aerosol density profiles in arbitrary units from ceilometer measurements between November 2009 and February 2013 are interpolated at measurement day time t and at retrieval grid middle heights h to provide a relative aerosol profile. This profile, turned into a partial optical depth per layer by multiplying with the layer thickness, is scaled to match the total aerosol extinction from AERONET τaer. The profile is then converted back into an intensive3 quantity by division by layer thickness, usually known as the aerosol extinction profile (extinction per unit length, AE). This is then used to calculate the aerosol scattering coefficient σaer in each layer by multiplying the layer aerosol extinction with the aerosol single-scattering albedo ωaer.

The aerosol absorption coefficient aaer is the layer aerosol extinction times (1−ωaer), and the aerosol phase function coefficients βaer,l can be calculated via the Henyey–Greenstein phase function and the asymmetry parameter g (see Eq. 5) e.g., , where l denotes the moment,

$\begin{array}{}\text{(5)}& {\mathit{\beta }}_{\mathrm{aer},l}=\left(\mathrm{2}l+\mathrm{1}\right)\cdot {g}^{l}.\end{array}$

### 4.2.3 Calculating final VLIDORT input

As mentioned in the introduction to this section, the layer input parameters for VLIDORT are total optical depth τ and single-scattering albedo ω in the layer, as well as the phase function coefficients βl. The quantities which have been calculated so far are the aerosol and air contributions to βl and the air and aerosol scattering coefficients σair and σaer, as well as the absorption coefficients agas from the trace gas in question and aaer from the aerosol in each VLIDORT input layer (red boxes in Fig. 3). These quantities need to be combined to total layer input parameters.

The total layer optical depth is simply the product of the layer thickness and the sum of all extinction and scattering coefficients; see Eq. (6):

$\begin{array}{}\text{(6)}& \mathit{\tau }=\mathrm{\Delta }h\cdot \left({\mathit{\sigma }}_{\mathrm{air}}+{\mathit{\sigma }}_{\mathrm{aer}}+{a}_{\mathrm{gas}}+{a}_{\mathrm{aer}}\right).\end{array}$

The combined single-scattering albedo can be calculated as in Eq. (7):

$\begin{array}{}\text{(7)}& \mathit{\omega }=\left({\mathit{\sigma }}_{\mathrm{air}}+{\mathit{\sigma }}_{\mathrm{aer}}\right)\cdot \mathrm{\Delta }h/\mathit{\tau }.\end{array}$

The combined expansion coefficients are calculated as in Eq. (8):

$\begin{array}{}\text{(8)}& {\mathit{\beta }}_{\text{l}}=\left({\mathit{\beta }}_{\mathrm{air},\mathrm{l}}\cdot {\mathit{\sigma }}_{\mathrm{air}}+{\mathit{\beta }}_{\mathrm{aer},\mathrm{l}}\cdot {\mathit{\sigma }}_{\mathrm{aer}}\right)/\left({\mathit{\sigma }}_{\mathrm{air}}+{\mathit{\sigma }}_{\mathrm{aer}}\right).\end{array}$

Since Jacobians with respect to changes in quantity χ in each layer are required, the normalized derivatives of the total primary optical quantities, τ, ω and β, with respect to this quantity need to be supplied. In the case of trace gas retrieval, χ=agasΔh, and in the case of aerosol retrieval, $\mathit{\chi }={\mathit{\sigma }}_{\mathrm{aer}}\mathrm{\Delta }h+{a}_{\mathrm{aer}}\mathrm{\Delta }h$. Therefore, what needs to be done is to calculate the normalized derivatives of Eqs. (6)–(8) with respect to these quantities. It should be remembered that β is a vector.

The three quantities for aerosol are calculated according to Eqs. (9)–(11).

$\begin{array}{ll}\text{(9)}& \frac{\mathit{\chi }}{\mathit{\tau }}\frac{\text{d}\mathit{\tau }}{\text{d}\mathit{\chi }}=& \frac{\left({\mathit{\sigma }}_{\mathrm{aer}}+{a}_{\mathrm{aer}}\right)\mathrm{\Delta }h}{\mathit{\tau }}\frac{\text{d}\mathit{\tau }}{\text{d}\left({\mathit{\sigma }}_{\mathrm{aer}}+{a}_{\mathrm{aer}}\right)\mathrm{\Delta }h}=& \frac{{\mathit{\sigma }}_{\mathrm{aer}}+{a}_{\mathrm{aer}}}{e}\\ \text{(10)}& \frac{\mathit{\chi }}{\mathit{\omega }}\frac{\text{d}\mathit{\omega }}{\text{d}\mathit{\chi }}=& \frac{\left({\mathit{\sigma }}_{\mathrm{aer}}+{a}_{\mathrm{aer}}\right)\mathrm{\Delta }h}{\mathit{\omega }}\cdot \left[\frac{\partial \mathit{\omega }}{\partial {\mathit{\sigma }}_{\mathrm{aer}}}\right& \frac{\text{d}{\mathit{\sigma }}_{\mathrm{aer}}}{\text{d}\left[\left({\mathit{\sigma }}_{\mathrm{aer}}+{a}_{\mathrm{aer}}\right)\mathrm{\Delta }h\right]}+\frac{\partial \mathit{\omega }}{\partial {a}_{\mathrm{aer}}}\frac{\text{d}{a}_{\mathrm{aer}}}{\text{d}\left[\left({\mathit{\sigma }}_{\mathrm{aer}}+{a}_{\mathrm{aer}}\right)\mathrm{\Delta }h\right]}]\\ & =\frac{{\mathit{\sigma }}_{\mathrm{aer}}{a}_{\mathrm{gas}}-{a}_{\mathrm{aer}}{\mathit{\sigma }}_{\mathrm{air}}}{\mathit{\omega }{e}^{\mathrm{2}}}\\ \text{(11)}& \frac{\mathit{\chi }}{\mathbit{\beta }}\frac{\text{d}\mathbit{\beta }}{\text{d}\mathit{\chi }}=& \frac{\left({\mathit{\sigma }}_{\mathrm{aer}}+{a}_{\mathrm{aer}}\right)\mathrm{\Delta }h}{\mathbit{\beta }}\frac{\text{d}\mathbit{\beta }}{\text{d}\left({\mathit{\sigma }}_{\mathrm{aer}}+{a}_{\mathrm{aer}}\right)\mathrm{\Delta }h}=& \frac{{\mathit{\sigma }}_{\mathrm{aer}}}{\mathbit{\beta }}\frac{{\mathbit{\beta }}_{\mathrm{aer}}-\mathbit{\beta }}{{\mathit{\sigma }}_{\mathrm{aer}}+{\mathit{\sigma }}_{\mathrm{air}}}\end{array}$

Here, e is the total extinction coefficient: $e=\mathit{\tau }/\mathrm{\Delta }h$

For the trace gas Jacobian calculation, the corresponding normalized derivatives are according to Eqs. (12)–(14).

$\begin{array}{}\text{(12)}& & \frac{\mathit{\chi }}{\mathit{\tau }}\frac{\mathrm{d}\mathit{\tau }}{\mathrm{d}\mathit{\chi }}=\frac{{a}_{\mathrm{gas}}\mathrm{\Delta }h}{\mathit{\tau }}\frac{\mathrm{d}\mathit{\tau }}{\mathrm{d}{a}_{\mathrm{gas}}}=\frac{{a}_{\mathrm{gas}}}{e}\text{(13)}& & \frac{\mathit{\chi }}{\mathit{\omega }}\frac{\mathrm{d}\mathit{\omega }}{\mathrm{d}\mathit{\chi }}=\frac{{a}_{\mathrm{gas}}\mathrm{\Delta }h}{\mathit{\omega }}\frac{\mathit{\omega }}{{a}_{\mathrm{gas}}}=-\frac{{a}_{\mathrm{gas}}}{e}\text{(14)}& & \frac{\mathit{\chi }}{\mathit{\beta }}\frac{\mathrm{d}\mathit{\beta }}{\mathrm{d}\mathit{\chi }}=\mathrm{0}\end{array}$

## 4.3 Calculating dSCDs and weighting functions

The forward model outputs, intensities and intensity Jacobians need to be converted into differential slant column densities and corresponding Jacobians. For each set of dSCDs, we have to run the forward model two times: Once with the gas absorption included and once without the gas absorption. Each of these sets consists of simulations in the desired telescope directions and one simulation towards the zenith. Note that the output from VLIDORT is Kx; hence the immediate output from VLIDORT has to be divided by x in case of linear retrieval to give the Jacobian alone.

In the following, the intensity and Jacobian of the simulation towards the desired angles without gas absorption are denoted ${I}_{\mathrm{0}}^{\mathit{\alpha }}$ and ${\mathbit{K}}_{\mathrm{0}}^{\mathit{\alpha }}$, the intensity (Jacobian) with gas absorption towards the zenith ${I}_{\mathrm{g}}^{\text{z}}$ (${\mathbit{K}}_{\mathrm{g}}^{\text{z}}$), the intensity (Jacobian) with gas absorption towards the desired angle ${I}_{\mathrm{g}}^{\mathit{\alpha }}$ (${\mathbit{K}}_{\mathrm{g}}^{\mathit{\alpha }}$) and the intensity (Jacobian) without gas absorption towards the zenith ${I}_{\mathrm{0}}^{\text{z}}$ (${\mathbit{K}}_{\mathrm{0}}^{\text{z}}$). The dSCD can be calculated as in Eq. (15).

$\begin{array}{}\text{(15)}& {\mathrm{dSCD}}^{\mathit{\alpha }}=\mathrm{ln}\left(\frac{{I}_{\mathrm{0}}^{\mathit{\alpha }}\cdot {I}_{\mathrm{g}}^{\text{z}}}{{I}_{\mathrm{g}}^{\mathit{\alpha }}\cdot {I}_{\mathrm{0}}^{\text{z}}}\right)/{\mathit{\xi }}_{\mathit{\lambda }}\end{array}$

The corresponding weighting function K is calculated as in Eq.16:

$\begin{array}{}\text{(16)}& {\mathbit{K}}^{\mathit{\alpha }}=\frac{\left({\mathbit{K}}_{\mathrm{0}}^{\mathit{\alpha }}{I}_{\mathrm{g}}^{\text{z}}{I}_{\mathrm{g}}^{\mathit{\alpha }}{I}_{\mathrm{0}}^{\text{z}}+{I}_{\mathrm{0}}^{\mathit{\alpha }}{\mathbit{K}}_{\mathrm{g}}^{\text{z}}{I}_{\mathrm{g}}^{\mathit{\alpha }}{I}_{\mathrm{0}}^{\text{z}}-{I}_{\mathrm{0}}^{\mathit{\alpha }}{I}_{\mathrm{g}}^{\text{z}}{\mathbit{K}}_{\mathrm{g}}^{\mathit{\alpha }}{I}_{\mathrm{0}}^{\text{z}}-{I}_{\mathrm{0}}^{\mathit{\alpha }}{I}_{\mathrm{g}}^{\text{z}}{I}_{\mathrm{g}}^{\mathit{\alpha }}{\mathbit{K}}_{\mathrm{0}}^{\mathit{\alpha }}\right)}{\left({I}_{\mathrm{0}}^{\mathit{\alpha }}{I}_{\mathrm{g}}^{\text{z}}{I}_{\mathrm{g}}^{\mathit{\alpha }}{I}_{\mathrm{0}}^{\text{z}}{\mathit{\xi }}_{\mathit{\lambda }}\right)}.\end{array}$
5 Error analysis

The errors of gas profile retrievals in remote-sensing applications can be divided into the following four different types (Rodgers2000):

1. smoothing error, arising from the constraint of the NO2 and aerosol profiles;

2. measurement noise error, arising from the noise in the spectra;

3. parameter error, including different parameters as aerosol profiles and cross sections; and the

4. forward model error, which originates from the simplification and uncertainties in the radiative transfer algorithm.

The gain matrix, as defined by Eq. (17), describes the change in retrieved x when measurement y changes. It can hence be used to map the uncertainty in measurement space to a state-vector uncertainty.

$\begin{array}{}\text{(17)}& \mathbf{\text{gain}}:=\frac{\partial \mathbit{x}}{\partial \mathbit{y}}={\left({\mathbf{K}}_{i}^{\text{T}}{\mathbf{S}}_{\mathrm{m}}^{-\mathrm{1}}{\mathbf{K}}_{i}+\mathbf{R}\right)}^{-\mathrm{1}}{\mathbf{K}}_{i}^{\text{T}}{\mathbf{S}}_{\mathrm{m}}\end{array}$

The error analysis in general produces an error pattern ϵy in the measurement space of vectors of the dSCDs retrieved by QDOAS and simulated by MMF. With help of the gain matrix, this error pattern is then mapped as ϵx into the space of the solution state.

## 5.1 Smoothing error

The retrieved atmospheric state vector x only represents a smoothed version of the true real state. The partial column averaging kernel matrix AKpcol describes how a change in the true atmospheric state vector xtrue is translated into changes in the retrieved state vector x, according to Eq. (18):

$\begin{array}{}\text{(18)}& \mathbit{x}-{\mathbit{x}}_{\text{a}}={\mathbf{AK}}_{\text{pcol}}\left({\mathbit{x}}_{\text{true}}-{\mathbit{x}}_{\text{a}}\right)+{\mathbit{ϵ}}_{x}.\end{array}$

The averaging kernel of the partial column and total column (AKtot) describe how the retrieved solution profile depends on the real atmosphere xtrue. They have to be taken into account if the profile and the column are used. The part which is not explained by AKpcol is explained by the retrieval error ϵx.

As MMF uses constrained least-square fitting, the AKpcol is calculated analytically by MMF itself for each scan according to Eq. (19):

The averaging kernel matrix is in general asymmetrical. The columns represent response functions related to a perturbation in a certain layer. The rows of the matrix represent the sensitivity of the NO2 partial column of a certain layer (or concentration in a certain layer) to all the different true partial columns or concentrations.

If AKpcol is expressed in units of partial columns (as here indicated by the subscript “pcol”), the sum of the rows is the total column averaging kernel. It represents the sensitivity of the vertical column density (VCD) to the anomalies in different heights.

The change in averaging kernel under transformation of the units is expressed by two matrix multiplications: one by multiplying the diagonal matrix containing the partial air column in its diagonal Uaircols from the right with the averaging kernel and the other by multiplying the inverse of the diagonal matrix ${\mathbf{U}}_{\text{aircols}}^{-\mathrm{1}}$ from the left side with the averaging kernel; see Eq. (20):

$\begin{array}{}\text{(20)}& {\mathbf{AK}}_{\text{VMR}}={\mathbf{U}}_{\text{aircols}}^{-\mathrm{1}}\cdot {\mathbf{AK}}_{\text{pcol}}\cdot {\mathbf{U}}_{\text{aircols}}.\end{array}$

The trace of the AKtot matrix represents the DOF, the number of independent pieces of information in the profile retrieval, which is around two for the NO2 retrieval.

Equation (18) separates the smoothing effect described with AK matrices from the other errors. Before comparing the retrieved profiles of MMF with profiles from models or other measurements, such as satellite measurements, it is necessary to smooth the profile from the other source with the AK from MMF . Alternatively, a smoothing error Ssmooth for the profile can be calculated according to from the a priori covariance matrix Sa and the averaging kernel. The latter is describing how the true atmospheric state varies, as a best estimate. Equation (21) describes how sensitivity of the retrieved atmospheric state vector depends on the true atmospheric state:

$\begin{array}{}\text{(21)}& \mathbf{AK}=\frac{\partial \mathbit{x}}{\partial {\mathbit{x}}_{\text{true}}}.\end{array}$

Both quantities can be given in different representations adjusted to the atmospheric state vector as a fraction of the a priori, the VMR (“VMR”), as number densities or as partial column profiles (“pcol”) as used above. Using AKVMR, Ssmooth can be calculated according to Eq. (22):

$\begin{array}{}\text{(22)}& {\mathbf{S}}_{\text{smooth}}=\left({\mathbf{AK}}_{\text{VMR}}-\mathbf{1}\right){\mathbf{S}}_{\text{a}},\text{VMR}\left({\mathbf{AK}}_{\text{VMR}}-\mathbf{1}{\right)}^{\text{T}}.\end{array}$

The smoothing error for the VCD is then calculated from the full covariance matrix of the smoothing error Ssmooth using a total column operator; see Eq. (23):

$\begin{array}{}\text{(23)}& {\mathit{\sigma }}_{\text{smooth}}=\sqrt{{\mathbit{g}}^{\text{T}}\cdot \underbrace{\left({\mathbf{AK}}_{\text{pcol}}-\mathbf{1}\right){\mathbf{S}}_{\text{a,pcol}}\left({\mathbf{AK}}_{\text{pcol}}-\mathbf{1}{\right)}^{\text{T}}{\mathbf{S}}_{\text{smooth}}^{\text{pcol}}\cdot \mathbit{g}}},\end{array}$

where g is the total column operator for partial column profiles: ${\mathbit{g}}^{\text{T}}=\left(\mathrm{1},\mathrm{1},\mathrm{1},\mathrm{1},\mathrm{1},\mathrm{1},\mathrm{1},\mathrm{1},\mathrm{1},\mathrm{\dots }\mathrm{1}\right)$.

As described before, the NO2 retrieval code uses a single a priori profile and covariance matrix taken from the chemical transport model WRF-Chem. This results in the sensitivities and smoothing errors represented in Fig. 4a.

The aerosol retrieval, however, uses a Tikhonov constraint and no covariance matrix is available. The a priori values for the aerosol extinction profile are reconstructed by the actual total aerosol optical depth reported by the daily AERONET measurements and the average vertical distribution, reconstructed from ceilometer measurements for each hour of the day. The covariance matrix for the aerosols Sa, aerosol is obtained by assuming a 100 % variability of the used a priori profile and an exponential correlation length of η=500 m between the different layers, as in the recent work by ; see Eq. (24):

$\begin{array}{}\text{(24)}& {\mathbf{S}}_{\text{a, aerosol}}\left[i,j\right]=\left(\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathit{%}\cdot \mathbit{A}{\mathbit{E}}_{\text{a}}\left[i\right]\cdot \mathbit{A}{\mathbit{E}}_{\text{a}}\left[j\right]\right){\mathrm{exp}}^{\left(-|z\left[i\right]-z\left[j\right]|/\mathit{\eta }\right)}.\end{array}$

Figure 4Averaging kernel columns: (a) in VMR representation for NO2 profile retrieval and (b) in the corresponding intensive quantity for the aerosol extinction (AE) profile. Each colored line shows the expected response on perturbation in a certain layer and belongs to the lower axis. The total column averaging kernel, (a) for the NO2 VCD and (b) for the total optical depth, is given as a black line and belongs to the upper axis. The vertical dashed line in each panel indicates an ideal sensitivity of 1.0 at all altitudes.

In the cases where no vertical aerosol extinction profile (AE profile) retrieval is available because the retrieval failed, or no O4 dSCDs were available in the first place, a NO2 profile can still be retrieved by using the a priori AE profile instead of the retrieved one. Since in that case no information of the O4 absorption is used, the estimated smoothing error Ssmooth is then directly given by the a priori covariance matrix Sa,aerosol for the AE profile. The a priori information about the optical properties described by the aerosol extinction profile is designed for cloud-free days and therefore the error analysis is just valid for such cloud-free days.

We have assumed a constant sensitivity AK and a constant a priori covariance matrix Sa for the calculation. The AK of the trace gas profile indeed depends strongly on the aerosol profile. It also depends slightly on the trace gas profile itself. The Sa covariance matrix of the aerosol extinction profile should be given for each hour as the a priori profile varies with time of day. Therefore, the estimation of the smoothing error as it is calculated here gives just a rough general idea about the size of the smoothing error.

## 5.2 Measurement noise

The measurement noise error can be calculated from the gain matrix and the measurement noise matrix (Sm) (Rodgers2000). The diagonal matrix (Sm) is already used in the optimal estimation of the profile retrieval to weight the different dSCDs of a scan according to the square of the errors from the QDOAS retrievals. The statistics of the measurement dSCD errors and its dependence on the elevation angle have been presented in Fig. 2.

The measurement noise error matrix for both (1) the NO2 profile and (2) aerosol extinction profile are calculated directly during the retrieval and are available from the MMF output of each measurement sequence as a full covariance matrix Snoise, given in units of partial columns and partial optical depths. The NO2 total column error is then calculated from this and the total column operator, and amounts to around 2.4 % of the VCD (see line 2 in Table 1). The errors in the profile are shown in Fig. 6.

Table 1Errors in NO2 vertical column density as a fraction of a retrieved VCD of 3.2×1016 molec. cm−2 measured at 13:58 on 20 May 2016 (UTC  6). The total error is calculated as the square root of the sum of squares of different independent errors. Different total errors are calculated, by including or not including the algorithm error , considering the two cases: that an aerosol extinction profile is retrieved from O4 absorption or just an a priori guess is used. For the assumed variability of the aerosol extinction, the error due to the uncertainty remaining after an aerosol extinction profile retrieval is 5.1 % instead of 9.8 %. However, if the algorithm error according to is included, the remaining error due to the uncertainty in the aerosol profile is slightly better: 9.3 % instead of the 9.8 % without O4 retrieval.

## 5.3 Parameter errors

The parameter errors originate from all uncertainties in input parameters in the forward model that are not properly fitted. These are, for example, the cross sections and, for the NO2 profile retrieval, also the used aerosol extinction profile.

### 5.3.1 Parameter error from spectroscopy

The error originating from an uncertainty in the cross section of 3 % is also around 3.0 % in the vertical column (see line 3 in Table 1) and similar in the vertical profile. This error can be calculated using the generalized measurement error covariance matrix ${\mathbf{S}}_{\text{spectroscopic}}^{y}$ in the measurement space assuming that all retrieved dSCDs are 3 % too high (or low) and using the outer product of the measurement state vector containing the retrieved dSCDs; see Eq. (25):

$\begin{array}{}\text{(25)}& {\mathbf{S}}_{\text{spectroscopic}}^{\mathbit{y}}=\left(\mathrm{3.0}\phantom{\rule{0.125em}{0ex}}\mathit{%}{\right)}^{\mathrm{2}}\cdot \mathbit{y}{\mathbit{y}}^{\text{T}}.\end{array}$

Using the gain matrix (see Eq. 17), this can be translated into an error in the retrieval space, as demonstrated in Eq. (26):

$\begin{array}{}\text{(26)}& {\mathbf{S}}_{\text{spectroscopic}}^{\mathbit{x}}=\mathbf{gain}\cdot {\mathbf{S}}_{\text{spectroscopic}}^{\mathbit{y}}\cdot {\mathbf{gain}}^{\text{T}}.\end{array}$

From the ${\mathbf{S}}_{\text{spectroscopic}}^{\mathbit{x}}$ matrix, the error in the total column and profile is calculated as shown earlier in Sect. 5.1. As expected, the error for the total column is 3 % and the errors in the NO2 VMR profile are shown in Fig. 6. The spectroscopic error is a purely systematic error and affects all retrievals in the same manner.

### 5.3.2 Parameter error from aerosol profile

The AE profile is crucial for the NO2 profile retrieval because of its strong contribution to the air mass factor. The uncertainties in this vector of input parameters arise from errors in the aerosol extinction profile retrieved from the spectral signature of the O4 dimer. For the propagation of the error of the O4 AE retrieval into the NO2 retrieval, we assume here that the total error in the AE profile depends on just the two contributions: (a) smoothing error and (b) measurement noise errors.

The effect of the assumed aerosol content in each layer on the retrieved NO2 profile is calculated by a sensitivity study. For this, 22 NO2 retrievals are performed for the same measurement sequence (dSCDs) using slightly modified aerosol extinction profiles. First, a normal retrieval with the best estimated aerosol extinction profile is retrieved. Then, the aerosol extinction profile is modified so that the partial extinction of the ith layer is increased by 1 % of the total optical depth with respect to the original aerosol extinction profile.

Figure 5Plots showing the results of how sensitive the NO2 profile is to changes in the aerosol extinction profile. (a) Columns of the matrix D created to describe the response in the NO2 VMR profile (the differences with respect to the original retrieval). (b) Same as (a) but the response to the perturbation in the aerosol profile is given as a fraction of the retrieved NO2 profile.

The difference between the perturbed and original NO2 VMR profiles is combined into the matrix D, thus describing how the NO2 VMR profile responds to changes in the aerosol at different heights. The result is presented in Fig. 5.

With the help of this matrix, DVMR, the different errors in the aerosol extinction profile, as measurement noise (see Sect. 5.2), smoothing error (see Sect. 5.1) or even due to the algorithm error according to , can be propagated to the NO2 profile and its VCD. The uncertainty in aerosol, represented by its Sa, aerosol matrix (see Eq. 24), translates into an uncertainty in NO2 VMR according to Eq. (27):

$\begin{array}{}\text{(27)}& {\mathbf{S}}_{{\mathrm{NO}}_{\mathrm{2}}-\text{VMR}}^{\text{aerosol}}={\mathbf{D}}_{\text{VMR}}{\mathbf{S}}_{\text{a,aerosol}}{\mathbf{D}}_{\text{VMR}}^{\text{T}}.\end{array}$

Using the technique from Eq. (23) with the unit transformation for the AK as in Eq. (20), this translates into an error on the VCD as in Eq. (28).

$\begin{array}{ll}\text{(28)}& & {\text{Error}}_{{\mathrm{NO}}_{\mathrm{2}}-\text{VCD}}^{\text{aerosol}}=& \sqrt{{g}^{\text{T}}{\mathbf{U}}_{\text{aircols}}\cdot {\mathbf{D}}_{\text{VMR}}\cdot {\mathbf{S}}_{\text{a,aerosol}}\cdot {{\mathbf{D}}_{\text{VMR}}}^{\text{T}}{\mathbf{U}}_{\text{aircols}}\cdot g}\end{array}$

The propagation of the smoothing and measurement noise errors of the O4 retrieval into the NO2 retrieval (separately resulting in 4.6 % and 2.2 %; see lines 6 and 5 in Table 1, respectively) results in a 5.1 % error in the NO2 VCD (see line 9 in Table 1). If the algorithm error introduced by (7.8 %, line 8 in Table 1) is included in the error propagated to NO2, the error in NO2 VCD is 9.3 % (line 10 in Table 1). However the algorithm error is calculated from the resulting residual of the fit and is dependent on the other error sources as mentioned earlier.

If no O4 retrieval is performed successfully, the aerosol a priori profile is used for the NO2 retrieval. The error in this aerosol profile is estimated by its aforementioned assumed variability and Sa covariance matrix: a variability of 100 % of the a priori profile is assumed for calculating the diagonal elements and an exponential correlation length of 500 m is used for calculating the off-axis elements of the Sa matrix. The error in the NO2 VCD would in this case be 9.8 % (see line 7 in Table 1).

Figure 6Altitude-dependent errors in a NO2 VMR profile measured at 13:58 on 20 May 2016 (UTC  6). (a) The square root of the diagonal in the covariance matrix describing the individual error contributions in the VMR profile and (b) the retrieved NO2 VMR profile with corresponding error bars (total error and variability). Panel (a) also includes the a priori profile in green.

## 5.4 Forward model error

The forward model error could be evaluated if an improved forward model was available. However, this is not the case. Errors in the forward model would result in systematic structures in the residual and a larger residual than expected from the error calculated by QDOAS for the slant columns. The least supported simplification in the forward model is the assumption of a horizontal homogeneity of gas and aerosols. A horizontal inhomogeneity leads to a set of slant columns in a scan which cannot be simulated by VLIDORT and the residuals of measured minus calculated slant columns indicate this error. Following , we calculate the variance in the residuals as a function of viewing angle and hence the algorithm error. But as all errors increase the residual the so-called algorithm error is not identical to the forward model error. However, it is maybe a good way to check empirically if there are important error contributions missing in the analytical analysis.

## 5.5 Total error

The results from the error estimations are summarized in Table 1. The overall error in the NO2 VCD is estimated to be around 14.1 % (20.3 % including the algorithm errors for NO2 and O4 calculated from the residuals). The contributions are 12.5 % from smoothing, 2.4 % from measurement noise, 3 % from spectroscopy and 5.1 % from errors in the aerosol extinction profile. The algorithm errors are 12.3 % from the NO2 retrieval itself and 7.8 % from the O4 AE retrieval. These errors are similar to those reported by .

The error in the vertical column is smaller than the errors in the VMR profile for almost all layers (Fig. 6). This can be explained by an anticorrelation in different partial column errors indicated by the full error covariance matrix. While the measurement noise seems to play a minor role, the smoothing and the aerosol profile are the main sources of error. Even though the error might seem large, the retrievals still provide new and relevant information of NO2 within the boundary layer of Mexico City.

6 Results

In this section, we present results of the NO2 variability measured at one of the stations comprising the MAX-DOAS network operated in Mexico City, and compare them with in situ measurements performed at the surface. The location of the MAX-DOAS station (see for a map) coincides with the AERONET station and the location of the surface sensor (19.32 N, 99.18 W). The launch location of the radiosondes is close to the airport, around 15 km northeast from UNAM at 19.40 N, 99.20 W.

A data set from January 2015 to July 2016 (approx. 19 months) is considered in this study, in which both MAX-DOAS and in situ surface measurements were available. A total number of 2531 coincidences of hourly averages were found in this period, considering that at least six measurement sequences of the MAX-DOAS instrument in each reported hour needed to be available. The complete time series is presented in Fig. 7, showing hourly NO2 averages of the in situ surface concentration in red, the total vertical column from the MAX-DOAS in blue and the average VMR of its first six layers closest to the surface in green. This plot includes all MAX-DOAS results regardless of the sky conditions.

Figure 7Hourly NO2 averages of in situ surface concentrations (red) in parts per billion, total VCDs from the MAX-DOAS instrument (blue) in molecules per cubic centimeter and the average VMR of the first six layers (up to 550 m) in the retrieved MAX-DOAS profiles (green), also in parts per billion.

Apart from the large data gaps in the beginning of 2016, the time series is quite complete and the fitted annual periodic functions (a Fourier series constrained to a fixed seasonal cycle shape) applied to the three data sets show a similar pattern. VMR values from the surface measurements are clearly higher than the VMRs detected in the lowest layers of the MAX-DOAS retrievals.

The 2531 coincidences were correlated averaging different numbers of layers starting from the ground and the resulting Pearson's coefficients (R) are plotted in Fig. 8a. Currently, the integration times in the spectra from which the O4 dSCDs are calculated are not long enough to ensure an O4 dSCD error resulting in DOF larger than 1 for the aerosol retrieval. Since we use a Tikhonov regularization for aerosol retrieval, this means that we can basically retrieve the total aerosol extinction. Due to the topographical circumstances in Mexico City, where the boundary layer rises steadily during the day, using an a priori profile that is calculated from hourly averages of ceilometer data will generally provide a good a priori profile shape. However, for days that likely have a very different aerosol profile from our a priori profile, the aerosol profile included in the trace gas forward model accounts for far more than the estimated error in Sect. 5.3.2 since this estimation uses a fixed percentage of the a priori profile in the Sa matrix. For cloudy days, for example, the form of the profile will look considerably different than the ceilometer averages. We are currently not able to retrieve such profiles. Hence, we limit our analysis to cloud-free days. Currently, the retrieval chain has no cloud-screening algorithm included. As an ad hoc solution, we use the presence of AERONET data on that day as an estimator for cloud-free days. For our entire data set, for about half of the coincident measurements (1270) there was at least one AERONET measurement available on that particular day while for the rest of our data (1261), no AERONET data were available on that day.

Figure 8Pearson's correlation coefficients (a) obtained from hourly surface and MAX-DOAS coincident measurements as a function of the number of layers considered in the VMR averages. The slopes obtained from the linear regressions, performed forcing the intercept to be at zero to the same data sets, are shown in (b). The layers below layer 6 are very likely inside the mixed boundary layer, while layers above layer 6 might be on the edge or above the boundary layer during some hours in the morning . The blue dashed lines only show data when AERONET data are available within 2 h of the scan and hence limit the measurements to cloud-free conditions. This selection criterion improves the correlation significantly.

Generally, a better correlation is found when averaged MAX-DOAS VMRs are compared to in situ surface measurements. Average R values larger than 0.7 are obtained in all cases when the first six to eight layers are considered in the correlation. These correlations decrease rapidly when fewer layers are considered due to increased errors in the lowermost part of the profile. Naturally, the R value also decreases towards a larger number of layers averaged since different air masses are measured at higher altitudes. As expected, slightly larger correlations are obtained for MAX-DOAS retrievals using a measured aerosol input from AERONET (blue traces in Fig. 8), which corresponds to clear-sky conditions.

Linear regressions such as the ones shown in Fig. 9 were generated for each data set, in order to gain more information as to how the number of layers averaged relates to the VMR comparison with the in situ surface measurements. The slope of all coincident data using the first seven to eight layers is 0.45, and it can reach values above 0.6 considering the first six layers in clear-sky conditions (AERONET data available within ±1 h of the measurement time). We thus decide to use the six-layer VMR averages of our MAX-DOAS retrievals in this study, reaching a height of 550 m above the ground level; in the comparisons we present with surface concentrations.

Figure 9Correlation plot of the MAX-DOAS six-layer VMR averages vs. the surface measurements shown separately for clear-sky measurements in which aerosol data on that day were available (blue) and for MAX-DOAS retrievals with no aerosol optical depth availability from AERONET, corresponding most likely to cloudy days (red).

To investigate the difference in our data set between clear and cloudy days, a correlation plot is presented in Fig. 9 showing the linear regressions produced from the MAX-DOAS six-layer VMR averages for the retrievals with AERONET availability on that day in blue and those retrievals using a default aerosol a priori profile in red. A significant improvement in the correlation coefficient going from 0.54 to 0.74 is evident in the plot. When all the coincident data are considered, regardless of if the retrieval had data available from the AERONET instrument on that day or not, the R and slope values are 0.62 and 0.39, respectively.

Examples of how the averaged six-layer VMRs and VCD compare to the ground level concentrations on individual days are presented in Fig. 10. These examples were chosen to depict distinct diurnal patterns that occur in different times of the year. Again, the MAX-DOAS six-layer product is consistently lower than the NO2 concentrations measured at the surface. These, together with the corresponding VCDs plotted on a different y axis, follow the pattern of the surface measurements quite well. The MAX-DOAS instrument captures the features observed at the surface but not without some interesting differences.

Figure 10Example days showing the NO2 variability in both VMR from the MAX-DOAS (green, <550 m) and surface measurements (red), as well as the total vertical column densities (blue).

Figure 11Diurnal NO2 variability from hourly means of the entire data set spanning from January 2015 to July 2016. Only the data with coincident surface and MAX-DOAS measurements are included for consistency in the comparison. Vertical bars correspond to the standard deviations.

Figure 12Annual NO2 variability resulting from the monthly mean coincident data (surface and MAX-DOAS) between January 2015 and July 2016. The standard deviations are given in vertical bars.

The data of 4 May 2016 was in the middle of a high-pollution episode in Mexico City, in which a restriction on the use of private vehicles was declared for 4 d. Ozone had surpassed, at least in one of the stations from the monitoring network operated by the city government, the 165 ppb 1 h average concentration. It is interesting to see in the 4 May plot that indeed the surface concentrations of nitrogen dioxide were high in the morning, a condition favoring ozone production, but there is a large NO2 peak just after noon in the MAX-DOAS VCD and six-layer product data which is not detected at the surface. This is more evident from the individual measurements shown as dots than from the hour averages.

Other differences have to do with the evolution of the mixed layer height, which has been shown to have a rapid growth in the late morning, strongly influencing the surface concentrations . This is for example evident in the 22 December plot, where a peak is observed at 11 h and a second one at 14 h. This peak is strongest at 11 h in the case of the surface measurement, but the total column shows the peak at 14 h to be dominant. The mixed layer has grown in the mean time so that the registered 14 h surface concentration is lower despite there being more NO2 molecules in the atmosphere.

In order to see the diurnal and seasonal patterns from both MAX-DOAS and surface measurements, all coincident measurements between January 2015 and July 2016 including days with clouds were averaged to a specific hour or month. For consistency, only the coincident 1 h data were considered in order to make the data sets comparable. The mean diurnal evolution with standard deviations as vertical bars is shown in Fig. 11. The maxima of both surface and MAX-DOAS six-layer products are at around 10 h, while that of the vertical column is shifted towards noon. This is expected from what is known from the emissions, growing mixed layer and ventilation patterns in Mexico City (e.g., García-Franco et al.2018). Despite the fact that the offset in the curves for surface- and MAX-DOAS measurements appears to be nearly constant throughout the day, it would be interesting to investigate further how this offset varies in different seasons, particularly when vertical mixing is not favored.

Figure 12, however, presents the compiled data considered in this study as monthly means in order to observe the seasonal variability and compare the data sets from these two measurement techniques. The three data products coincide in having the highest monthly mean values between April and June. A difference which is probably interesting to note is during December and January, when the concentrations from the surface measurements are relatively higher than those from the MAX-DOAS products. This might have to do with the fact that during these winter months the mixed layer is shallower than during the rest of the year, and the sensitivity of the MAX-DOAS instrument is reduced for the lowest layers as described above.

7 Conclusions

In this contribution we describe the methodology used to analyze the data produced by the network of MAX-DOAS instruments operating in the Mexico City metropolitan area. In general terms, MMF is a retrieval code developed to process the data acquired by the instruments converting the measured spectra together with other input parameters to vertical profiles of aerosols and certain trace gases. We have tested the performance of the code for NO2 at one of the stations and present some diagnostics and a full error analysis of the results. In the case of NO2, all error sources amount to about 14.1 %, the smoothing error being the dominating one (12.5 %) followed by the error due to the uncertainties in the aerosol extinction profile (5.1 %).

Both the resulting total vertical column densities and a product consisting of the average VMR in the first six layers (<550 m above the ground level) were compared to ground level NO2 concentrations. Good correlations were obtained between the six-layer product and the values from the surface measurements, with R values between 0.6 and 0.7. However, the MAX-DOAS systematically underestimates the ground level concentrations by a factor of about 0.6. This is consistent with the total column averaging kernels and variability Sa reported in Sect. 5.1, meaning that the MAX-DOAS instrument has a significantly lower sensitivity near the surface and is most sensitive at a height of around 1 km. It is shown, however, that this underestimation is less for clear-sky conditions as suggested from comparing data when aerosol optical depths were available from independent measurements.

Although results are shown only from one instrument located in the southern part of the city (UNAM), several years of data are available from three other stations at different locations within the metropolitan area (Acatlán, Vallejo and Cuautitlán) that are being analyzed and used for studying the spatial and temporal variability of NO2 in conjunction with several satellite products.

In this work we present a new and competitive retrieval code which has been developed for retrieving trace gas profiles from MAX-DOAS measurements. Since this study, it has been further improved to perform retrievals in logarithmic space to avoid unphysical negative partial columns and oscillations; it uses the more stable Levenberg–Marquardt iteration scheme and decouples the retrieval and simulation grids4.

Further efforts include the improvement in the measurement noise by increasing the integration times during spectral acquisition. This will allow us to improve the retrieval of other gases such as formaldehyde (HCHO) and other weak absorbers. It will also enable us to use OE as a retrieval strategy for aerosol retrieval and possibly increase the DOF to comparable values as for the NO2 retrieval. This, together with an inclusion of a cloud-screening algorithm in the retrieval chain, will make the retrieval less dependent on the availability of AERONET data.

Data availability
Data availability.

The MAX-DOAS data from the UNAM station can be accessed via http://www.epr.atmosfera.unam.mx/maxdoas_data/no2/nidforval/ (last access: 9 April 2019, ). Data from other stations, with other versions or periods measured within Mexico City's MAX-DOAS network, should be requested from Claudia Rivera (claudia.rivera@atmosfera.unam.mx). Measurements from the surface monitor (CCA station) are performed by Mexico City's Automatic Monitoring Network and can be accessed via http://www.aire.cdmx.gob.mx/ (last access: 9 April 2019).

Appendix A

In a recent update of the code, implemented after the analysis presented here (i.e., not used for obtaining the results here), the retrieval space was changed from linear space to logarithmic retrieval space. This means that the retrieval works in a linear (dSCD measurement) and logarithmic (profile retrieval) space now. This enhances the nonlinearity of the problem and requires a change of iteration scheme. The GN iteration scheme (Eq. 1) is replaced by a slightly slower but more stable Levenberg–Marquardt (LM) iteration scheme; i.e., Eq. (1) is replaced by Eq. (A1) for more nonlinear inversion problems (Rodgers2000).

$\begin{array}{ll}\text{(A1)}& {\mathbit{x}}_{i+\mathrm{1}}=& \phantom{\rule{0.25em}{0ex}}{\mathbit{x}}_{i}+{\left[\left(\mathrm{1}+\mathit{\gamma }\right)\mathbf{R}+{\mathbf{K}}_{i}^{\text{T}}\right]{\mathbf{S}}_{\mathrm{m}}^{-\mathrm{1}}{\mathbf{K}}_{i}\right]}^{-\mathrm{1}}& \left[{\mathbf{K}}_{i}^{\text{T}}{\mathbf{S}}_{\mathrm{m}}^{-\mathrm{1}}\left(\mathbit{y}-\mathbit{F}\left({\mathbit{x}}_{i}\right)\right)-\mathbf{R}\left({\mathbit{x}}_{i}-{\mathbit{x}}_{\text{a}}\right)\right]\end{array}$

The symbols mean the same as in Eq. (1): superscript “T” denotes transposed, superscript “−1” denotes the inverse and subscript “a” indicates a priori values. Sm is the measurement error covariance matrix, y denotes the vector of measured differential slant column densities, F(xi) denotes the simulated differential slant column densities at iteration i, calculated using the forward model with input profile xi, and Ki is the Jacobian matrix at the ith iteration.

The new xi+1 in Eq. (A1) is only accepted if the cost function in Eq. (A2) decreases with respect to the cost function of the previous iteration.

$\begin{array}{ll}\text{(A2)}& \text{cost}=& \phantom{\rule{0.25em}{0ex}}\sum _{k}^{\text{angles}}\sum _{j}^{\text{angles}}\left(\mathbit{y}-\mathbit{F}\left(\mathbit{x}\right){\right)}_{k}{{\mathbf{S}}_{\mathrm{m}}^{-\mathrm{1}}}^{kj}\left(\mathbit{y}-\mathbit{F}\left(\mathbit{x}\right){\right)}_{j}& +\sum _{k}^{\text{layers}}\sum _{j}^{\text{layers}}\left(\mathbit{x}-{\mathbit{x}}_{\text{a}}{\right)}_{k}\mathbf{S}{{a}^{-\mathrm{1}}}^{kj}\left(\mathbit{x}-{\mathbit{x}}_{\text{a}}{\right)}_{j}\end{array}$

If this is the case, xi+1 is accepted and if (1+γ) is not yet equal to 1, the factor (1+γ) is halved for the next iteration. If, however, the cost function increases, the newly calculated xi+1 is discarded and the ith calculation repeated with a factor (1+γ) increased by a factor of 16.

In order to counteract the slowdown of the retrieval, more restrictions were placed on the observation geometry for a single scan: a single relative azimuth angle and a single solar zenith angle per scan. This means in particular that two different viewing directions cannot be treated as a single scan any longer. Also, it is no longer necessary to have the full simulation grid as the retrieval grid. The retrieval can be limited to the lower part of the simulation grid. Although the former change (with respect to observation geometry) means a significant cut in flexibility, it results in a retrieval time speed up of a factor of 4, and a more typical retrieval time per scan and per constituent (i.e., for both aerosol and trace gas) is around 5 s.

Tests using the logarithm of the partial layer vertical column density (for NO2 retrieval) or layer extinction profile (for aerosol retrieval) motivated the change to the LM iteration scheme due to the increased nonlinearity when working in a semi-log space as state-measurement space. With this new configuration, MMF has been participating in the round robin comparison of different retrieval codes for the FRM4DOAS project . It has also participated in the profile retrieval from dSCD from the CINDI-2 campaign, both for NO2 and HCHO .

If the retrieval is to be performed in log space, i.e., to work with ln (x) instead of with x as the retrieval parameter, K in Eq. (16) is multiplied by x (since the raw output from VLIDORT is Kx, no division by x is necessary because the quantity needed in logarithmic space is Kx) and the uncertainty covariance matrix Sa in Eq. (A1) needs to be left and right multiplied by 1∕x.

The LM scheme of Eq. (A1) has currently only been tested with OE and not with Tikhonov regularization; i.e., the aerosol retrieval was also performed using OE.

Appendix B: Calculations for dry air density

In this appendix, we outline the calculation of Z, the molar mass of air Ma and the depolarization ratio Δ. We follow . If T is given in kelvin (temperature in degrees Celsius is denoted as TC), σ the wavenumber per micrometer, pressure P in pascals and volume mixing ratio of CO2 in parts per million, the constants in Eqs. (B1)–(B14) have the following values.

$\begin{array}{llr}{k}_{\mathrm{0}}& =& \mathrm{238.0185}\\ {a}_{\mathrm{0}}& =& \mathrm{1.58123}×{\mathrm{10}}^{-\mathrm{6}}\\ {b}_{\mathrm{0}}& =& \mathrm{5.707}×{\mathrm{10}}^{-\mathrm{6}}\\ {w}_{\mathrm{0}}& =& \mathrm{295.235}\\ {k}_{\mathrm{1}}& =& \mathrm{5792105.0}\\ {a}_{\mathrm{1}}& =& -\mathrm{2.9331}×{\mathrm{10}}^{-\mathrm{8}}\\ {b}_{\mathrm{1}}& =& -\mathrm{2.051}×{\mathrm{10}}^{-\mathrm{8}}\\ {w}_{\mathrm{1}}& =& \mathrm{2.6422}\\ {k}_{\mathrm{2}}& =& \mathrm{57.362}\\ {a}_{\mathrm{2}}& =& \mathrm{1.1043}×{\mathrm{10}}^{-\mathrm{10}}\\ {c}_{\mathrm{0}}& =& \mathrm{1.9898}×{\mathrm{10}}^{-\mathrm{4}}\\ {w}_{\mathrm{2}}& =& -\mathrm{0.03238}\\ {k}_{\mathrm{3}}& =& \mathrm{167917.0}\\ \mathit{\alpha }& =& \mathrm{1.00062}\\ {c}_{\mathrm{1}}& =& -\mathrm{2.376}×{\mathrm{10}}^{-\mathrm{6}}\\ {w}_{\mathrm{3}}& =& \mathrm{0.004028}\\ d& =& \mathrm{1.83}×{\mathrm{10}}^{-\mathrm{11}}\\ \mathit{\beta }& =& \mathrm{3.14}×{\mathrm{10}}^{-\mathrm{8}}\\ A& =& \mathrm{1.2379}×{\mathrm{10}}^{-\mathrm{5}}\\ C& =& \mathrm{33.937}\\ e& =& -\mathrm{0.765}×{\mathrm{10}}^{-\mathrm{8}}\\ \mathit{\gamma }& =& \mathrm{5.6}×{\mathrm{10}}^{-\mathrm{7}}\\ B& =& -\mathrm{1.912}×{\mathrm{10}}^{-\mathrm{2}}\\ D& =& -\mathrm{6.343}×{\mathrm{10}}^{\mathrm{3}}\end{array}$

First, the saturation vapor pressure 𝒮 is calculated according to Eq. (B1), then the enhancement factor for water vapor in air, , is calculated according to Eq. (B2) and with this, the molar fraction of water vapor in moist air xw is calculated according to Eq. (B3), where Pw is partial water vapor pressure in pascals.

$\begin{array}{}\text{(B1)}& & \mathcal{S}=\mathrm{exp}\left(A\phantom{\rule{0.125em}{0ex}}{T}^{\mathrm{2}}+B\phantom{\rule{0.125em}{0ex}}T+C+D/T\right)\text{(B2)}& & \mathcal{F}=\mathit{\alpha }+\mathit{\beta }\phantom{\rule{0.125em}{0ex}}P+\mathit{\gamma }\phantom{\rule{0.125em}{0ex}}{T}_{\text{C}}^{\mathrm{2}}\text{(B3)}& & {\mathbit{x}}_{\text{w}}=\mathcal{F}\phantom{\rule{0.125em}{0ex}}{P}_{\text{w}}/P\end{array}$

The air density in dry air ${n}_{\mathrm{air},\mathrm{dry},\mathrm{450}}$ (0 % humidity) at standard conditions (15 C, 101.325 Pa, xCO2=450 ppm) can be calculated as in Eq. (B4) and from this the dry air density at xCO2 is calculated according to Eq. (B5).

$\begin{array}{}\text{(B4)}& \left({n}_{\mathrm{air},\mathrm{dry},\mathrm{450}}-\mathrm{1}\right)×{\mathrm{10}}^{\mathrm{8}}=\frac{{k}_{\mathrm{1}}}{{k}_{\mathrm{0}}-{\mathit{\sigma }}^{\mathrm{2}}}+\frac{{k}_{\mathrm{3}}}{{k}_{\mathrm{2}}-{\mathit{\sigma }}^{\mathrm{2}}}\end{array}$

$\begin{array}{ll}\text{(B5)}& \left({n}_{\mathrm{air},\mathrm{dry},{x\mathrm{CO}}_{\mathrm{2}}}-\mathrm{1}\right)=& \phantom{\rule{0.25em}{0ex}}\left({n}_{\mathrm{air},\mathrm{dry},\mathrm{450}}-\mathrm{1}\right)×{\mathrm{10}}^{\mathrm{8}}& \left(\mathrm{1.0}+\mathrm{0.534}×{\mathrm{10}}^{-\mathrm{6}}\cdot \left(x{\mathrm{CO}}_{\mathrm{2}}-\mathrm{450.0}\right)\right)\end{array}$

The equivalent quantity for water vapor, nwv, is calculated as follows:

$\begin{array}{}\text{(B6)}& {n}_{\mathrm{wv}}=\mathrm{1.022}\cdot \left({w}_{\mathrm{0}}+{w}_{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathit{\sigma }}^{\mathrm{2}}+{w}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathit{\sigma }}^{\mathrm{4}}+{w}_{\mathrm{3}}{\mathit{\sigma }}^{\mathrm{6}}\right).\end{array}$

The molar mass of dry air, Ma, with xCO2 is calculated as

$\begin{array}{}\text{(B7)}& {M}_{\text{a}}={\mathrm{10}}^{-\mathrm{3}}\left(\mathrm{28.9635}+\mathrm{12.011}×{\mathrm{10}}^{-\mathrm{6}}\left(x{\mathrm{CO}}_{\mathrm{2}}-\mathrm{400.0}\right)\right).\end{array}$

The compressibility of dry air, Za, and pure water vapor, Zw, are calculated as

$\begin{array}{ll}\text{(B8)}& {Z}_{\text{a}}=& \phantom{\rule{0.25em}{0ex}}\mathrm{1.0}-\frac{\mathrm{101325.0}}{\mathrm{288.15}}\phantom{\rule{0.125em}{0ex}}\left({a}_{\mathrm{0}}+\mathrm{15}\phantom{\rule{0.125em}{0ex}}{a}_{\mathrm{1}}+\mathrm{225}\phantom{\rule{0.125em}{0ex}}{a}_{\mathrm{2}}\right),\text{(B9)}& {Z}_{\text{w}}=& \phantom{\rule{0.25em}{0ex}}\mathrm{1.0}-\frac{\mathrm{1333.0}}{\mathrm{293.15}}\phantom{\rule{0.125em}{0ex}}\left({a}_{\mathrm{0}}+\mathrm{20}\phantom{\rule{0.125em}{0ex}}{a}_{\mathrm{1}}+\mathrm{400}\phantom{\rule{0.125em}{0ex}}{a}_{\mathrm{2}}+{b}_{\mathrm{0}}\right& +\mathrm{20}\phantom{\rule{0.125em}{0ex}}{b}_{\mathrm{1}}+{c}_{\mathrm{0}}+\mathrm{20}\phantom{\rule{0.125em}{0ex}}{c}_{\mathrm{1}})+{\left(\frac{\mathrm{1333.0}}{\mathrm{293.15}}\right)}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\left(d+e\right).\end{array}$

Using Eq. (B10) once with the values for Z, R and T for standard water vapor and once with the values for standard air, one can calculate ρws and ρaxs, respectively.

$\begin{array}{}\text{(B10)}& \mathit{\rho }=\left(p{M}_{\text{a}}/\left(Z{R}_{\mathrm{gas}}T\right)\left(\mathrm{1}-{\mathbit{x}}_{\text{w}}\left(\mathrm{1}-{M}_{\mathrm{mwv}}/{M}_{\text{a}}\right)\right)\end{array}$

Here, Rgas=8.3144621 J mol−1 K−1 is the gas constant and Mmwv=0.018015 kg mol−1 is the molar mass of water vapor. Z, the compressibility of moist air at T, xw and P can be calculated as

$\begin{array}{ll}\text{(B11)}& Z=& \phantom{\rule{0.25em}{0ex}}\mathrm{1.0}-\left(P/T\right)\left({a}_{\mathrm{0}}+{a}_{\mathrm{1}}{T}_{\text{C}}+{a}_{\mathrm{2}}{T}_{\text{C}}^{\mathrm{2}}+\left({b}_{\mathrm{0}}+{b}_{\mathrm{1}}{T}_{\text{C}}\right){\mathbit{x}}_{\text{w}}& +\left({c}_{\mathrm{0}}+{c}_{\mathrm{1}}{T}_{\text{C}}\right){\mathbit{x}}_{\text{w}}^{\mathrm{2}}\right)+\left(P/T{\right)}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\left(d+e{\mathbit{x}}_{\text{w}}^{\mathrm{2}}\right).\end{array}$

Using again Eq. (B10), with the actual values for P, xw and Z, the density of the dry air component can be calculated according to Eq. (B12) and the density of the water vapor component can be calculated according to Eq. (B13).

$\begin{array}{}\text{(B12)}& & {\mathit{\rho }}_{\mathrm{dryair}}=P{M}_{\text{a}}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{1}-{\mathbit{x}}_{\text{w}}\right)/\left(Z\phantom{\rule{0.125em}{0ex}}{R}_{\mathrm{gas}}\phantom{\rule{0.125em}{0ex}}T\right)\text{(B13)}& & {\mathit{\rho }}_{\mathrm{wv}}=P\phantom{\rule{0.125em}{0ex}}{M}_{\mathrm{mwv}}{\mathbit{x}}_{\text{w}}/\left(Z\phantom{\rule{0.125em}{0ex}}{R}_{\mathrm{gas}}\phantom{\rule{0.125em}{0ex}}T\right)\end{array}$

Finally, the number density of the air ρ and the refractive index n are calculated as

$\begin{array}{}\text{(B14)}& & \mathit{\rho }=\frac{{\mathit{\rho }}_{\mathrm{dryair}}}{{M}_{\text{a}}}\mathit{Av},\text{(B15)}& & n=\frac{{\mathit{\rho }}_{\mathrm{dryair}}}{{\mathit{\rho }}_{\mathrm{axs}}}\left({n}_{\text{air,dry,}x{\mathrm{CO}}_{\mathrm{2}}}-\mathrm{1}\right)+\frac{{\mathit{\rho }}_{\mathrm{wv}}}{{\mathit{\rho }}_{\mathrm{ws}}}\left({n}_{\mathrm{wv}}-\mathrm{1}\right).\end{array}$

Here, $\mathit{Av}=\mathrm{6.02214129}×{\mathrm{10}}^{\mathrm{23}}$ is the Avogadro number.

We use the N2, O2 and CO2 depolarization factors from and fixed VMR as summarized in Table B1 and calculate the wavelength-dependent depolarization factor via Eq. (B16), where F is the so-called King factor and is calculated for N2 (λ in micrometers) and O2 according to Eqs. (B18) and (B19), respectively. The complete depolarization factor is then calculated according to Eq. (B19).

Table B1Depolarization ratios and VMR for main air constituents

$\begin{array}{ll}\text{(B16)}& & \mathrm{\Delta }=\mathrm{6.0}\cdot \left(\mathrm{1}-F\right)/\left(-\mathrm{7}F-\mathrm{3}\right)\text{(B17)}& & F\left({\mathrm{N}}_{\mathrm{2}}\right)=\mathrm{1.034}+\mathrm{3.17}×{\mathrm{10}}^{-\mathrm{4}}/{\mathit{\lambda }}^{\mathrm{2}}\text{(B18)}& & F\left({\mathrm{O}}_{\mathrm{2}}\right)=\mathrm{1.096}+\mathrm{1.385}×{\mathrm{10}}^{-\mathrm{3}}/{\mathit{\lambda }}^{\mathrm{2}}& +\mathrm{1.448}×{\mathrm{10}}^{-\mathrm{4}}/{\mathit{\lambda }}^{\mathrm{4}}\\ \text{(B19)}& & \mathrm{\Delta }=\sum \left({\mathrm{\Delta }}_{i}{\mathrm{VMR}}_{i}\right),\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}i={\mathrm{N}}_{\mathrm{2}},{\mathrm{O}}_{\mathrm{2}},\mathrm{Ar},{\mathrm{CO}}_{\mathrm{2}}\end{array}$

For the Rayleigh scattering cross section, we implement Eq. (7.37) from and Eq. (19.5) from ; see Eq. (B20).

$\begin{array}{}\text{(B20)}& {Q}_{\mathrm{Ray}}=\frac{\mathrm{24}{\mathit{\pi }}^{\mathrm{3}}\left({n}^{\mathrm{2}}-\mathrm{1}{\right)}^{\mathrm{2}}}{{\mathit{\lambda }}^{\mathrm{4}}{\mathit{\rho }}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\left({n}^{\mathrm{2}}+\mathrm{2}{\right)}^{\mathrm{2}}}\frac{\mathrm{6}+\mathrm{3}\mathrm{\Delta }}{\mathrm{6}-\mathrm{7}\phantom{\rule{0.125em}{0ex}}\mathrm{\Delta }}\end{array}$
Author contributions
Author contributions.

MMF is responsible for the retrieval code development and testing, the retrieval chain setup from spectra to profiles, and the retrieval parameter choices for MMF and software support. In terms of text contributions, MMF wrote Sect. 4, Appendix A and B and part of the abstract and Sect. 1. MMF produced Figs. 1–3. CR is responsible for the QDOAS retrieval setup and parameter choices. CR wrote Sect. 3 and parts of Sect. 2. WS and ZO are responsible for the data analysis and interpretation, as well as for the sensitivity tests and the error analysis. They wrote Sect. 5 and parts of Sects. 6 and 7. They produced Figs. 4–12. JA constructed and installed the instruments and was responsible for their maintenance. JA wrote Sect. 2. AB provided technical support for instruments and data management. JAGR provided the WRF-Chem simulations that serve as a priori data for NO2 retrieval. JAGR wrote parts of Sect. 4.2.1. MG was involved in the data analysis and interpretation and wrote Sect. 7 and parts of Sects. 1 and 6.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

We would like to thank the financial support provided by CONACYT (grants 275239 and 290589) and DGAPA-UNAM (grants IN107417, IN111418, and IA100716 and Martina M. Friedrich's postdoctoral fellowship). SEDEMA is acknowledged for providing the in situ surface concentrations from their monitoring network (RAMA) as well for partially supporting the MAX-DOAS network, and Arne Krueger, Alfredo Rodríguez, Hector Soto and Delibes Flores are thanked for their technical assistance. We thank Caroline Fayt (caroline.fayt@aeronomie.be), Michel Van Roozendael (michelv@aeronomie.be) and Thomas Danckaert for the free use of the QDOAS software and we thank Robert Spurr for free use of the VLIDORT radiative transfer code package. We thank the Mexican Solarimetric Service for their effort in establishing and maintaining the AERONET Mexico City site. We thank the University of Wyoming Department of Atmospheric Science for providing the sounding data on http://weather.uwyo.edu/upperair/sounding.html (last access: 9 April 2019).

This research has been supported by the CONACYT (grant no. 275239), the CONACYT (grant no. 290589), the DGAPA-UNAM (grant no. IN107417), the DGAPA-UNAM (grant no. IN111418) and the DGAPA-UNAM (grant no. IA100716).

Review statement
Review statement.

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

References

Arellano, J., Krüger, A., Rivera, C., Stremme, W., Friedrich, M., Bezanilla, A., and Grutter, M.: The MAX-DOAS network in Mexico City to measure atmospheric pollutants, Atmosfera, 29, 157–167, 2016. a, b, c, d, e

Bates, D. R.: Rayleigh scattering by air, Planet. Space Sci., 32, 785–790, https://doi.org/10.1016/0032-0633(84)90102-8, 1984. a, b

Buchwitz, M., Rozanov, V. V., and Burrows, J. P.: Development of a correlated-k distribution band model scheme for the radiative transfer program GOMETRAN/SCIATRAN for retrieval of atmospheric constituents from SCIMACHY/ENVISAT-1 data, Proc. SPIE, 3495, Satellite Remote Sensing of Clouds and the Atmosphere III, https://doi.org/10.1117/12.332681, 1998. a

Burrows, J. P., Richter, A., Dehn, A., Deters, B., Himmelmann, S., Voigt, S., and Orphal, J.: Atmospheric remote-sensing reference data from GOME. 2. Temperature dependent absorption cross sections of O3 in the 231-794 nm range, J. Quant. Spectrosc. Ra., 61, 509–517, https://doi.org/10.1016/S0022-4073(98)00037-5, 1999. a, b

Carlson, D., Donohoue, D., Platt, U., and Simpson, W. R.: A low power automated MAX-DOAS instrument for the Arctic and other remote unmanned locations, Atmos. Meas. Tech., 3, 429–439, https://doi.org/10.5194/amt-3-429-2010, 2010. a

Ciddor, P. E.: Refractive index of air: new equations for the visible and near infrared, Appl. Opt., 35, 1566, https://doi.org/10.1364/AO.35.001566, 1996. a

Coburn, S., Dix, B., Sinreich, R., and Volkamer, R.: The CU ground MAX-DOAS instrument: characterization of RMS noise limitations and first measurements near Pensacola, FL of BrO, IO, and CHOCHO, Atmos. Meas. Tech., 4, 2421–2439, https://doi.org/10.5194/amt-4-2421-2011, 2011. a

Dahlback, A. and Stamnes, K.: A new spherical model for computing the radiation field available for photolysis and heating at twilight, Planet. Space Sci., 39, 671–683, https://doi.org/10.1016/0032-0633(91)90061-E, 1991. a

Danckaert, T., Fayt, C., Van Roozendael, M., De Smedt, I., Letocart, V., Merlaud, A., and Pinardi, G.: QDOAS Software user manual, Belgian Institute for Space Aeronomy, 2013. a, b, c, d, e

Finlayson-Pitts, B. J. and Pitts, J. N.: Chemistry of the upper and lower atmosphere: theory, experiments, and applications, Academic Press, San Diego, CA, 2000. a

Frieß, U., Monks, P. S., Remedios, J. J., Rozanov, A., Sinreich, R., Wagner, T., and Platt, U.: MAX-DOAS O4 measurements: A new technique to derive information on atmospheric aerosols: 2. Modeling studies, J. Geophys. Res.-Atmos., 111, D14203, https://doi.org/10.1029/2005JD006618, 2006. a

Frieß, U., Beirle, S., Alvarado Bonilla, L., Bösch, T., Friedrich, M. M., Hendrick, F., Piters, A., Richter, A., van Roozendael, M., Rozanov, V. V., Spinei, E., Tirpitz, J.-L., Vlemmix, T., Wagner, T., and Wang, Y.: Intercomparison of MAX-DOAS Vertical Profile Retrieval Algorithms: Studies using Synthetic Data, Atmos. Meas. Tech. Discuss., https://doi.org/10.5194/amt-2018-423, in review, 2018. a

Galle, B., Johansson, M., Rivera, C., Zhang, Y., Kihlman, M., Kern, C., Lehmann, T., Platt, U., Arellano, S., and Hidalgo, S.: Network for Observation of Volcanic and Atmospheric Change (NOVAC) A global network for volcanic gas monitoring: Network layout and instrument description, J. Geophys. Res., 115, D05304, https://doi.org/10.1029/2009JD011823, 2010. a

García-Franco, J., Stremme, W., Bezanilla, A., Ruiz-Angulo, A., and Grutter, M.: Variability of the Mixed-Layer Height Over Mexico City, Bound.-Lay. Meteorol., 167, 493–507, https://doi.org/10.1007/s10546-018-0334-x, 2018. a, b, c, d, e

Goody, R. M. and Yung, Y. L.: Atmospheric radiation: theoretical basis, Oxford University Press, 1989. a, b

Hendrick, F., Müller, J.-F., Clémer, K., Wang, P., De Mazière, M., Fayt, C., Gielen, C., Hermans, C., Ma, J. Z., Pinardi, G., Stavrakou, T., Vlemmix, T., and Van Roozendael, M.: Four years of ground-based MAX-DOAS observations of HONO and NO2 in the Beijing area, Atmos. Chem. Phys., 14, 765–781, https://doi.org/10.5194/acp-14-765-2014, 2014. a

Hermans, C., Vandaele, A. C., Carleer, M., Fally, S., Colin, R., Jenouvrier, A., Coquart, B., and Mérienne, M. F.: Absorption cross-sections of atmospheric constituents: NO2, O2, and H2O, Environ. Sci. Pollut. Res., 6, 151–158, 1999. a, b

Hess, M., Koepke, P., and Schult, I.: Optical Properties of Aerosols and Clouds: The Software Package OPAC., B. Am. Meteorol. Soc., 79, 831–844, https://doi.org/10.1175/1520-0477(1998)079<0831:OPOAAC>2.0.CO;2, 1998. a

Hönninger, G., Bobrowski, N., Palenque, E. R., Torrez, R., and Platt, U.: Reactive bromine and sulfur emissions at Salar de Uyuni, Bolivia, J. Geophys. Res. Lett., 31, L04101, https://doi.org/10.1029/2003GL018818, 2004. a, b

Johansson, M., Rivera, C., de Foy, B., Lei, W., Song, J., Zhang, Y., Galle, B., and Molina, L.: Mobile mini-DOAS measurement of the outflow of NO2 and HCHO from Mexico City, Atmos. Chem. Phys., 9, 5647–5653, https://doi.org/10.5194/acp-9-5647-2009, 2009. a

Kirmse, B., Delon, A., and Jost, R.: NO2 absorption cross section and its temperature dependence, J. Geophys. Res., 102, 16089–16098, https://doi.org/10.1029/97JD00075, 1997. a

Kurucz, R. L., Furenlid, I., Brault, J., and Testerman, L.: Solar flux atlas from 296 to 1300 nm, National Solar Observatory Atlas No. 1, 1984. a

Marquard, L.: Modellierung des Strahlungstransports in der Erdatmosphäre für absorptionsspektroskopische Messungen im ultravioletten und sichtbaren Spektralbereich, Logos Verlag, Berlin, ISBN 3897220520, 1998. a

Marquard, L. C., Wagner, T., and Platt, U.: Improved air mass factor concepts for scattered radiation differential optical absorption spectroscopy of atmospheric species, J. Geophys. Res., 105, 1315–1327, https://doi.org/10.1029/1999JD900340, 2000. a

Melamed, M. L., Basaldud, R., Steinbrecher, R., Emeis, S., Ruíz-Suárez, L. G., and Grutter, M.: Detection of pollution transport events southeast of Mexico City using ground-based visible spectroscopy measurements of nitrogen dioxide, Atmos. Chem. Phys., 9, 4827–4840, https://doi.org/10.5194/acp-9-4827-2009, 2009. a

Meller, R. and Moortgat, G. K.: Temperature dependence of the absorption cross sections of formaldehyde between 223 K and 323 K in the wavelength range 225–375 nm, J. Geophys. Res., 105, 7089–7101, 2000. a

Molina, L. T. and Molina, M. J.: Air quality in the Mexico megacity: an integrated assessment, Kluwer Academic Publishers, Dordrecht, Boston, 2002. a

Orphal, J.: A critical review of the absorption cross-sections of O3 and NO2 in the ultraviolet and visible, J. Photoch. Photobio. A, 157, 185–209, https://doi.org/10.1016/S1010-6030(03)00061-3, 2002. a

Ortega, I., Koenig, T., Sinreich, R., Thomson, D., and Volkamer, R.: The CU 2-D-MAX-DOAS instrument – Part 1: Retrieval of 3-D distributions of NO2 and azimuth-dependent OVOC ratios, Atmos. Meas. Tech., 8, 2371–2395, https://doi.org/10.5194/amt-8-2371-2015, 2015. a

Palazzi, E., Petritoli, A., Giovanelli, G., Kostadinov, I., Bortoli, D., Ravegnani, F., and Sackey, S. S.: PROMSAR: A backward Monte Carlo spherical RTM for the analysis of DOAS remote sensing measurements, Adv. Space Res., 36, 1007–1014, https://doi.org/10.1016/j.asr.2005.05.017, 2005. a

Perliski, L. M. and Solomon, S.: On the evaluation of air mass factors for atmospheric near-ultraviolet and visible absorption spectroscopy, J. Geophys. Res., 98, 10363, https://doi.org/10.1029/93JD00465, 1993. a

Platt, U. and Stutz, J.: Differential Optical Absorption Spectroscopy, Springer-Verlag Berlin Heidelberg, https://doi.org/10.1007/978-3-540-75776-4, 2008. a, b

Platt, U., Pfeilsticker, K., and Vollmer, M.: Radiation and Optics in the Atmosphere, Springer Handbook of Lasers and Optics, p. 1165, https://doi.org/10.1007/978-0-387-30420-5_19, 2007. a, b

Rahpoe, N., von Savigny, C., Weber, M., Rozanov, A. V., Bovensmann, H., and Burrows, J. P.: Error budget analysis of SCIAMACHY limb ozone profile retrievals using the SCIATRAN model, Atmos. Meas. Tech., 6, 2825–2837, https://doi.org/10.5194/amt-6-2825-2013, 2013a. a

Rahpoe, N., von Savigny, C., Weber, M., Rozanov, A. V., Bovensmann, H., and Burrows, J. P.: Error budget analysis of SCIAMACHY limb ozone profile retrievals using the SCIATRAN model, Atmos. Meas. Tech., 6, 2825–2837, https://doi.org/10.5194/amt-6-2825-2013, 2013b. a

Rivera, C., Barrera, H., Grutter, M., Zavala, M., Galle, B., Bei, N., Li, G., and Molina, L. T.: NO2 fluxes from Tijuana using a mobile mini-DOAS during Cal-Mex 2010, Atmos. Environ., 70, 532–539, https://doi.org/10.1016/j.atmosenv.2012.12.026, 2013. a

Rivera, C., Stremme, W., Ojeda, Z., Bezanilla, A., Arellano, J., Friedrich, M. M., and Grutter, M.: NO2 profiles from MAXDOAS measurements at UNAM station, available at: http://www.epr.atmosfera.unam.mx/maxdoas_data/no2/nidforval/, last access: 9 April 2019. a

Rodgers, C. D.: Inverse Methods for Atmospheric Sounding: Theory and Practice, Series on Atmospheric, Oceanic and Planetary Physics, Vol. 2, World Scientific Publishing Co, 256 pp., https://doi.org/10.1142/3171, 2000. a, b, c, d, e

Rodgers, C. D. and Connor, B. J.: Intercomparison of remote sounding instruments, J. Geophys. Res.-Atmos., 108, 4116, https://doi.org/10.1029/2002JD002299, 2003. a

Rozanov, V. V., Rozanov, A. V., Kokhanovsky, A. A., and Burrows, J. P.: Radiative transfer through terrestrial atmosphere and ocean: Software package SCIATRAN, J. Quant. Spectrosc. Ra., 133, 13–71, https://doi.org/10.1016/j.jqsrt.2013.07.004, 2014. a

Rozanov, V. V., Dinter, T., Rozanov, A. V., Wolanin, A., Bracher, A., and Burrows, J. P.: Radiative transfer modeling through terrestrial atmosphere and ocean accounting for inelastic processes: Software package SCIATRAN, J. Quant. Spectrosc. Ra., 194, 65–85, https://doi.org/10.1016/j.jqsrt.2017.03.009, 2017. a

Spurr, R.: VLIDORT: A linearized pseudo-spherical vector discrete ordinate radiative transfer code for forward model and retrieval studies in multilayer multiple scattering media, J. Quant. Spectrosc. Ra., 102, 316–342, 2006. a, b

Spurr, R.: User's Guide VLIDORT Version 2.6, RT Solutions, Inc., 2013. a, b, c

Spurr, R. J. D., Kurosu, T. P., and Chance, K. V.: A linearized discrete ordinate radiative transfer model for atmospheric remote-sensing retrieval, J. Quant. Spectrosc. Ra., 68, 689–735, https://doi.org/10.1016/S0022-4073(00)00055-8, 2001. a, b, c

Stamnes, K., Tsay, S.-C., Jayaweera, K., and Wiscombe, W.: Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media, Appl. Opt., 27, 2502–2509, https://doi.org/10.1364/AO.27.002502, 1988. a

Steck, T.: Methods for determining regularization for atmospheric retrieval problems, Appl. Opt., 41, 1788–1797, https://doi.org/10.1364/AO.41.001788, 2002. a

Tirpitz, J.-L., Frieß, U., Hendrick, F., Alberti, C., Allaart, M., Apituley, A., Bais, A., Beirle, S., Benavent, N., Berkhout, S., Bösch, T., Bognar, K., Bruchkouski, I., Chan, K. L., Chengxin, Z., den Hoed, M., Donner, S., Drosoglou, T., Friedrich, M. M., Frumau, A., Gast, L., Gomez, L., Gielen, C., Hao, N., Hensen, A., Henzing, B., Hoque, S., Irie, H., Jin, J., Koenig, T. K., Kreher, K., Kuhn, J., Kumar, V., Lampel, J., Li, A., Liu, C., Ma, J., Merlaud, A., Mishra, A., Nieto, D., Peters, E., Pinardi, G., Piters, A., Pöhler, D., Postylyakov, O., Richter, A., van Roozendael, M., Schmitt, S., Sinha, V., Spinei, E., Stein, D., Swart, D., Tack, F., Vlemmix, T., van der Hoff, R., Vonk, J., Wagner, T., Wang, S., Wang, Y., Wang, Z., Wenig, M., Wiegner, M., Wittrock, F., Xie, P., Xing, C., Xu, J., and Zhao, X.: DOAS Vertical Profile Retrieval Algorithms: Studies on Field Data from the CINDI-2 Campaign, Atmos. Meas. Tech., in preparation, 2019. a

Vandaele, A. C., Hermans, C., Simon, P. C., Carleer, M., Colin, R., Fally, S., Mérienne, M. F., Jenouvrier, A., and Coquart, B.: Measurements of the NO2 absorption cross-section from 42 000 cm−1 to 10 000 cm−1 (238–1000 nm) at 220 K and 294 K, J. Quant. Spectrosc. Ra., 59, 171–184, 1998. a, b

Vandaele, A. C., Hermans, C., Fally, S., Carleer, M., Merienne, M.-F., Jenouvrier, A., Coquart, B., and Colin, R.: Absorption cross-sections of NO2: simulation of temperature and pressure effects, J. Quant. Spectrosc. Ra., 76, 373–391, https://doi.org/10.1016/S0022-4073(02)00064-X, 2003. a

Vidot, J., Jourdan, O., Kokhanosvky, A. A., Szczap, F., Giraud, V., and Rozanov, V. V.: Retrieval of tropospheric NO2 columns from satellite measurements in presence of cirrus: A theoretical sensitivity study using SCIATRAN and prospect application for the A-Train, J. Quant. Spectrosc. Ra., 111, 586–601, 2010.  a

von Friedeburg, C.: Derivation of Trace Gas Information combining Differential Optical Absorption Spectroscopy with Radiative Transfer Modelling, PhD thesis, Ruperto-Carola University of Heidelberg, Germany, 2003. a

von Friedeburg, C., Wagner, T., Geyer, A., Kaiser, N., Vogel, B., Vogel, H., and Platt, U.: Derivation of tropospheric NO3 profiles using off-axis differential optical absorption spectroscopy measurements during sunrise and comparison with simulations, J. Geophys. Res.-Atmos., 107, 4168, https://doi.org/10.1029/2001JD000481, 2002. a

Wagner, T., Beirle, S., Brauers, T., Deutschmann, T., Frieß, U., Hak, C., Halla, J. D., Heue, K. P., Junkermann, W., Li, X., Platt, U., and Pundt-Gruber, I.: Inversion of tropospheric profiles of aerosol extinction and HCHO and NO2 mixing ratios from MAX-DOAS observations in Milano during the summer of 2003 and comparison with independent data sets, Atmos. Meas. Tech., 4, 2685–2715, https://doi.org/10.5194/amt-4-2685-2011, 2011. a

Wang, S., Cuevas, C. A., Frieß, U., and Saiz-Lopez, A.: MAX-DOAS retrieval of aerosol extinction properties in Madrid, Spain, Atmos. Meas. Tech., 9, 5089–5101, https://doi.org/10.5194/amt-9-5089-2016, 2016. a

Wang, Y., Lampel, J., Xie, P., Beirle, S., Li, A., Wu, D., and Wagner, T.: Ground-based MAX-DOAS observations of tropospheric aerosols, NO2, SO2 and HCHO in Wuxi, China, from 2011 to 2014, Atmos. Chem. Phys., 17, 2189–2215, https://doi.org/10.5194/acp-17-2189-2017, 2017. a, b, c, d, e, f, g, h

Wilmouth, D. M., Hanisco, T. F., Donaue, N. M., and Anderson, J. G.: Fourier Transform Ultraviolet Spectroscopy of the ${A}^{\mathrm{2}}{\prod }_{\mathrm{3}/\mathrm{2}}←{X}^{\mathrm{2}}{\prod }_{\mathrm{3}/\mathrm{2}}$, J. Phys. Chem. A, 103, 8935–8945, 1999. a

Wittrock, F., Oetjen, H., Richter, A., Fietkau, S., Medeke, T., Rozanov, A., and Burrows, J. P.: MAX-DOAS measurements of atmospheric trace gases in Ny-Ålesund – Radiative transfer studies and their application, Atmos. Chem. Phys., 4, 955–966, https://doi.org/10.5194/acp-4-955-2004, 2004. a

In a recent update, the GN scheme was replaced by the more stable Levenberg–Marquardt (LM) iteration scheme; more details on recent changes can be found in Appendix A.

We only use the LIDORT part of VLIDORT, meaning that we only consider the total intensity and not the polarization.

Bulk property that does not change when changing the size of the system.

Details can be found in Appendix A.