BOREAS – a new MAX-DOAS profile retrieval algorithm for aerosols and trace gases

We present a new MAX-DOAS profiling algorithm for aerosols and trace gases, BOREAS, which utilizes an iterative solution method including Tikhonov regularization and the optimal estimation technique. The aerosol profile retrieval is based on a novel approach in which the absorption depth of O₄ is directly used in order to retrieve extinction coefficient profiles instead of the commonly used perturbation theory method. The retrieval of trace gases is done with the frequently used optimal estimation method but significant improvements are presented on how to deal with wrongly weighted a priori constraints and for scenarios in which the a priori profile is inaccurate.

Performance tests are separated into two parts. First, we address the general sensitivity of the retrieval to the example of synthetic data calculated with the radiative transfer model SCIATRAN. In the second part of the study, we demonstrate BOREAS profiling accuracy by validating the results with the help of ancillary measurements carried out during the CINDI-2 campaign in Cabauw, the Netherlands, in 2016.

The synthetic sensitivity tests indicate that the regularization between measurement and a priori constraints is insufficient when knowledge of the true state of the atmosphere is poor. We demonstrate a priori pre-scaling and extensive regularization tests as a tool for the optimization of retrieved profiles. The comparison of retrieval results with in situ, ceilometer, NO₂ lidar, sonde and long-path DOAS measurements during the CINDI-2 campaign always shows high correlations with coefficients greater than 0.75. The largest differences can be found in the morning hours, when the planetary boundary layer is not yet fully developed and the concentration of trace gases and aerosol, as a result of a low night-time boundary layer having formed, is focused in a shallow, near-surface layer.

1 Introduction

Aerosols and trace gases play an important role for life on Earth as high concentrations have adverse impacts on human health and climate. Furthermore, aerosols impact the Earth's energy budget by radiative forcing. They interact with solar radiation by scattering and absorption processes. Additionally, they have an impact on the formation of clouds when acting as cloud condensation nuclei (CCN). Despite considerable efforts by the scientific community, uncertainties in aerosol radiative forcing are still high (Stocker, 2014).

In addition to aerosols, atmospheric trace gases are of importance in an increasingly urbanized world as they impact human health, agriculture, acid deposition and climate (e.g. Lelieveld et al., 2002). Trace gases in the troposphere such as nitrogen dioxide and sulfur dioxide (NO₂, SO₂), nitrous acid (HONO), formaldehyde (HCHO) and glyoxal (CHOCHO) need to be measured to assess their impact on air quality and because their spatial and temporal inhomogeneous distributions provide information on different emission sources, e.g. the combustion of fossil fuels, biomass burning and agriculture (Hendrick et al., 2014; Irie et al., 2011; Pinardi et al., 2013; Wittrock et al., 2004). The separation and identification of anthropogenic from natural sources, as well as measurements of their temporal variability, improve the understanding of their effects on climate and life on Earth. The current lack of knowledge leads to uncertainties in temporal and spatial emission patterns of trace gases and aerosols and limits our ability to fully understand atmospheric processes. This results in a need for comprehensive measurements to fill these gaps.

For more than 15 years, Multi-Axis Differential Optical Absorption Spectroscopy (MAX-DOAS) measurements have been used to investigate the chemical composition of the troposphere (Bobrowski et al., 2003; Hönninger and Platt, 2002; Leser et al., 2003; Wagner et al., 2004; Wittrock et al., 2004; van Roozendael et al., 2004). This passive remote-sensing method is based on absorption spectroscopy applied to measurements of scattered sunlight. The advantages of ground-based MAX-DOAS when compared to satellite observations are the high sensitivity for the lowermost layers of the troposphere, the high temporal and spatial resolution of measurements and the lower cost.

The strength of the absorption signal detected in scattered sunlight depends on the absorber amounts and their vertical distribution but also on the length of the light path. In general, this length in a certain altitude layer is a function of the measurement geometry. Therefore, a set of MAX-DOAS radiance measurements taken at different elevation angles (lines of sight, LOSs) contains information about the vertical distribution of trace gases, which can be retrieved. However, the retrieval of absorber profiles from MAX-DOAS measurements is an ill-posed problem that needs additional constraints (see Sect. 3.2). Thus, the inversion, which is applied to retrieve vertical concentration profiles, is done by elaborated mathematical methods such as optimal estimation (OE). In addition to the trace gas profile, the aerosol extinction profile needs to be retrieved as well as it has a non-linear effect on the trace gas retrieval and it is too variable to be approximated by climatologies. In contrast to the trace gas retrievals, aerosol retrievals are strongly non-linear, necessitating iterative inversion schemes such as the Gauss–Newton algorithm or the Levenberg–Marquardt algorithm (Rodgers, 2004).

Profile retrieval algorithms used in the scientific community are either inversion algorithms or parameterized approaches. Inversion algorithms directly link the measurement quantity to the vertical profile of the target absorber with the help of a forward model (Clémer et al., 2010; Frieß et al., 2006; Hendrick et al., 2004; Wagner et al., 2004; Wang et al., 2017; Wittrock, 2006; Wittrock et al., 2004). This model is usually calculated with a radiative transfer model (RTM), assuming a selected set of atmospheric conditions. For trace gases, the model also includes a sensitivity matrix which can be considered derivative of the measurement quantity with respect to the trace gas concentration for each altitude level. For aerosols, this sensitivity matrix is normally calculated via perturbation theory. This perturbation approach is shortly summarized as follows: by consideration of a trace gas with a well-known vertical distribution such as the oxygen dimer O₂−O₂ (or short O₄), the aerosol extinction for each layer is changed gradually and the resulting modelled observations for different LOSs are compared with measurements until the iteration converges (e.g. the difference between the measurement and modelled quantity is small enough).

Parameterization, on the other hand, means that a forward model, based on a limited set of parameters, is used to describe the measurement quantity. Frequently used forward model parameters are integrated values and profile information such as shape and height of certain absorber layers (Lee et al., 2009; Li et al., 2010; Sinreich et al., 2005, 2013; Vlemmix et al., 2011; Wagner et al., 2011). The forward model results are calculated with an RTM and are least-squares fitted to the measurements. Generally, look-up tables (LUTs) are pre-calculated for a set of different scenarios, parameters and geometries covering all relevant atmospheric conditions, avoiding high computational efforts during near-real-time calculations.

Inversion algorithms have the advantage that they are not limited to the scenarios used when creating the LUT but unrealistic profiles are possible when the measurement, inversion or regularization (weighting between the information from measurements and a priori information) is poor. On the other hand, parameterized approaches evaluate profiles much faster, as the slow forward model computation was already done when creating the LUT.

Although efforts to derive concentration profiles from MAX-DOAS measurements have been made for more than one decade, profiling is still a difficult task and the results of different algorithms can differ strongly when the absorber of interest is highly variable and inhomogeneous on spatial and temporal scales (Frieß et al., 2018; Vlemmix et al., 2015; Zieger et al., 2011). Comparison studies of different profiling algorithms for synthetic as well as for real data, summarizing the current state of the art and including results from the algorithm introduced here, will be reported by Frieß et al. (2018) and Tirpitz et al. (2019).

The purpose of this study is the introduction of IUP Bremen's new MAX-DOAS profile retrieval algorithm for aerosols and trace gases BOREAS (Bremen Optimal estimation REtrieval for Aerosols and trace gaseS), which has been developed to improve the earlier profile retrieval algorithm (Wittrock, 2006). BOREAS uses a novel approach for the retrieval of aerosols but a similar optimal estimation technique for the retrieval of trace gases. In contrast to perturbation-based inversion algorithms, BOREAS uses the change (from an a priori state) of depth in an absorption band to get information on the aerosol content which caused this change. To the authors' knowledge, this approach has never been used within an operational profiling algorithm and it complements the variety of methods with another promising technique. In order to improve the standard OE method for trace gases, a simple way of changing the weighting between a priori and measurement constraints is introduced by varying a regularization factor so that oscillations are minimized. Additionally, the usage of a priori pre-scaling is highlighted as one way of improving the results when the a priori profile is far away from the real atmosphere.

The development of BOREAS is aimed at several key properties:

1. Flexibility: the algorithm should retrieve aerosol and trace gas profiles from any MAX-DOAS measurement with pre-filtering options for the data.

2. Accuracy and stability: the algorithm should be stable in terms of varying atmospheric conditions when retrieving profiles for several years of data. The profiling results (modelled observation) should fit the measured observations with a high accuracy.

3. Automation: the retrieval should respond to problematic data/settings (e.g. low information content, wrong regularization) automatically with an included problem solution scheme.

4. Fast: the algorithm should be fast enough to allow near-real-time profile retrievals.

The paper is organized as follows: in the first section, a typical MAX-DOAS measurement and its information content are shortly described. The next section focuses on the theoretical background of the retrieval algorithm. This section is followed by a detailed sensitivity study. The fourth section is divided into an analysis of synthetic data, a discussion of error sources and an example of application to real measurements, which shows that BOREAS is able to retrieve profiles with high reliability and accuracy. The final section (Sect. 5) concludes and summarizes the results presented in this study.

2 MAX-DOAS measurements and the DOAS retrieval

Modern ground-based MAX-DOAS instruments are capable of measuring scattered sunlight in the ultraviolet (UV) and visible (VIS) spectral range with a full azimuthal and elevation angle coverage of the hemisphere (see Fig. 1).

Figure 1 Schematic representation of typical MAX-DOAS measurement geometries.

In general, the spectral least-squares fit of absorption cross sections and a polynomial to the slant optical thickness (logarithm of a Sun-normalized radiance at a certain LOS, also called slant optical thickness, SOT) $\mathit{\tau }=-\ln (I/I_{\mathrm{0}})$ (DOAS fit) provides integrated concentrations ρ along the light path L, the slant column densities S (or SCD) (Platt and Stutz, 2008). The length of the light path L depends on the Sun's position, the scattering properties of the atmosphere and on the viewing geometry. For ground-based measurements, the extraterrestrial solar spectrum I₀ is usually replaced by an intensity spectrum measured in the zenith direction I_ref (reference spectrum) to compensate for Fraunhofer lines and instrumental issues (e.g. temperature dependent wavelength shifts). As this zenith spectrum itself is usually contaminated by the absorber of interest to some extent, the difference between the two SCDs, ΔS is retrieved and referred to as differential SCD. The differential optical thickness Δτ_j of a certain absorber j only, is defined as

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }{\mathit{\tau }}_j(\mathit{\lambda },\mathbf{\Omega })={\mathit{\sigma }}_j\left(\mathit{\lambda }\right)\mathrm{\Delta }{S}_j\left(\mathbf{\Omega }\right),\end{array}$$

where σ_j(λ) is the spectral absorption cross section and $\mathbf{\Omega }=\mathit{\{}\mathit{\alpha },\mathit{\phi },\mathit{\vartheta }\mathit{\}}$ represents the angular variables. Here, α is the elevation angle (or LOS), φ is the relative azimuth angle (RAA) between the Sun and viewing azimuth, and ϑ is the solar zenith angle (SZA).

When using a single zenith spectrum measured around noon as I_ref, the tropospheric signal in I_ref is minimized, but for absorbers, which are also present in the stratosphere, stratospheric absorption contributes significantly to ΔS in the morning and evening. This is due to diurnal variations in the stratospheric amounts of some absorbers such as NO₂ and the change in light paths through the stratosphere because of varying SZA (see e.g. Leser et al., 2003). To minimize the stratospheric signal, a reference spectrum measured close in time to the intensity I for a certain geometry is applied. Usually, a scan of measurements in different elevation angles is followed by a measurement in the zenith direction so that the whole scan can be fitted with the same reference spectrum.

The ratio of slant column density S and the vertically integrated absorber concentration $V={\int }_{\mathrm{0}}^H\mathit{\rho }\left(z\right)\text{d}z$ (vertical column density or VCD, V) is called the air mass factor, M (AMF).

$$\begin{array}{}\text{(2)}& M\left(\mathbf{\Omega }\right)\equiv {\displaystyle \frac{S\left(\mathbf{\Omega }\right)}{V}}\end{array}$$

Generally, this value depends on the radiative transfer through the atmosphere and on parameters such as SZA, surface albedo, wavelength and the profiles of aerosols and gaseous absorbers. For trace gases in the boundary layer, it can, in a first approximation, be described by the geometric air mass factor

$$\begin{array}{}\text{(3)}& M_{\text{geom}}\left(\mathit{\alpha }\right)={\displaystyle \frac{S\left(\mathit{\alpha }\right)}{V}}={\displaystyle \frac{\mathrm{1}}{\sin \left(\mathit{\alpha }\right)}},\end{array}$$

assuming that the trace gas layer is located below the last scattering point (Hönninger et al., 2004). The differential AMF (dAMF), assuming a zenith reference spectrum, is then $\mathrm{\Delta }M=S/V-S_{\text{ref}}/V=(\mathrm{1}-\sin (\mathit{\alpha }\left)\right)/\sin \left(\mathit{\alpha }\right)$. The difference between geometric and true AMF depends strongly on the wavelength, elevation angle, relative azimuth angle and aerosols.

Because of the altitude-dependent concentration of absorbers, there is a need for a similar quantity describing the ratio of S and V in a certain layer. The concept of the box air mass factor ℳ_z (BAMF) introduces this height dependence by defining the ratio for each altitude layer z:

$$\begin{array}{}\text{(4)}& {\mathcal{M}}_z\left(\mathbf{\Omega }\right)={\displaystyle \frac{S_z\left(\mathbf{\Omega }\right)}{V_z}}.\end{array}$$

Here, ℳ_z depends on the altitude in contrast to the total AMF M. Box air mass factors are calculated by radiative transfer models and can be understood as sensitivity of the partial slant column to the partial vertical column at a certain altitude. On the other hand, multiplication of a profile of partial vertical columns with ℳ_z leads to the corresponding slant column density.

3 Retrieval algorithm description

The profile retrieval algorithm BOREAS was developed in order to retrieve aerosol and trace gas vertical profiles from MAX-DOAS measurements. The aerosol retrieval is fully implemented within the RTM SCIATRAN (Rozanov et al., 2014) to decrease computation time for the iterative minimization scheme. BOREAS is a Python-written analysis script, which calls SCIATRAN for the aerosol retrieval and for calculations of BAMF matrices, which are then used within an optimal-estimation-based trace gas retrieval. A flow chart depicting BOREAS is shown in Fig. 2. In the next subsections, we give an overview of the individual steps of the algorithm.

Figure 2 BOREAS analysis flow chart. The algorithm is separated into two steps. Step 1: aerosol retrieval within a Tikhonov regularization. Step 2: optimal-estimation-based trace gas retrieval.

3.1 Retrieval of aerosol profiles

The standard DOAS fit does not provide direct information about aerosols present in the atmosphere. However, the scattering and absorption properties of aerosols have an impact on the measured differential slant column density ΔS(Ω) because scattering processes can significantly modify the light path. Since the single-scattering albedo ω (SSA) stays more or less constant when the aerosol type does not change over time and altitude, it is not a quantity to be retrieved in BOREAS and is kept constant with typical values for the prevalent aerosol type (e.g. urban pollution: SSA =0.92). The angular scattering distribution of aerosols is fully quantified by the scattering phase function. Within BOREAS, we use the Henyey–Greenstein approximation (Henyey and Greenstein, 1941) with constant values for the asymmetry factor g, which quantifies the amount of forward and backward scattering (e.g. urban pollution: g=0.68). In this parameterization, the optical properties of aerosols are fully defined with g, SSA and an extinction coefficient profile σ_e(z). The latter is the retrieval parameter for BOREAS aerosol retrievals. Note that the usage of measured phase functions will be implemented in BOREAS in the future, since Henyey–Greenstein is in some situations not an accurate representation of the atmospheric aerosol scattering distribution.

If σ_e(z) is known, then the aerosol optical thickness (AOT) can be determined as ${\mathit{\tau }}_{\text{aer}}={\int }_{\mathrm{0}}^H{\mathit{\sigma }}_{\text{e}}\left(z\right)\text{d}z$, where the integration is performed over the entire atmosphere. If the vertical profile of an absorber number density is known, differences between modelled and measured S for this absorber are the result of differences between the assumed and real aerosol profiles. As observations at different LOSs have varying sensitivity to the presence of aerosols at different altitudes, this can be used to retrieve an aerosol profile.

Generally, the vertical profiles of species in the troposphere are unknown because of the temporal and spatial variability of emission sources, transport and conversion processes. However, the oxygen monomer O₂ is in this respect an exceptional species because it only depends on pressure and temperature. Furthermore, as the oxygen dimer O₄ concentration is proportional to the squared monomer concentration, its profile exponentially decreases with altitude as well. The O₄ slant column density can easily be determined because O₄ has spectral absorption features in the wavelength regions of most DOAS fits. Detailed sensitivity studies of O₄ ΔS(Ω) measurements regarding changes in atmospheric aerosol properties can be found in Wagner et al. (2004) and Frieß et al. (2006).

In the BOREAS aerosol retrieval algorithm, the difference between modelled and measured O₄ differential slant optical thicknesses Δτ(λ,Ω) (DSOT) is used to retrieve aerosol extinction profiles in an iterative process. The measured Δτ(λ,Ω) is calculated using Eq. (1), where ΔS_i(Ω) is the DSCD retrieved in the framework of the standard DOAS fit, and the simulated O₄ DSOT ($\mathrm{\Delta }\stackrel{\mathrm{̃}}{\mathit{\tau }}(\mathit{\lambda },\mathbf{\Omega })$) is calculated as follows:

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }\stackrel{\mathrm{̃}}{\mathit{\tau }}(\mathit{\lambda },\mathbf{\Omega },N_{\text{a}}(z\left)\right)=\ln {\displaystyle \frac{I(\mathit{\lambda },\mathbf{\Omega },N_{\text{a}}(z\left)\right)}{I_{\text{ref}}(\mathit{\lambda },N_{\text{a}}(z\left)\right)}},\end{array}$$

where $I(\mathit{\lambda },\mathbf{\Omega },N_{\text{a}}(z\left)\right)$ and I_ref(λ,N_a(z)) are the intensities calculated using SCIATRAN under the assumption that O₄ is the only absorber. The dependency on the reference geometry Ω_ref is summarized with the index “ref” and will be neglected from now on. The vertical profiles of pressure and temperature are used according to the US standard atmosphere model (NASA, 1976) but can be replaced by measured atmospheric conditions when available. The aerosol loading is described using an a priori concentration number density vertical profile N_a(z).

The inverse problem with respect to the aerosol optical depth is then formulated as

$$\begin{array}{}\text{(6)}& {∥\mathrm{\Delta }\mathit{\tau }(\mathit{\lambda },\mathbf{\Omega })-\mathrm{\Delta }\stackrel{\mathrm{̃}}{\mathit{\tau }}(\mathit{\lambda },\mathbf{\Omega },N_{\text{a}}(z\left)\right)-P(\mathit{\lambda },\mathbf{\Omega })∥}^{\mathrm{2}}⟶min,\end{array}$$

where P(λ,Ω) is a polynomial of lower order and its argument Ω emphasizes that polynomial coefficients depend on the LOS angle. The assumptions made in this formulation are that (1) the O₄ absorption derived from the measurements does not depend on the concentration of other trace gases and (2) that the optical depth in an atmosphere without other absorbers can be described as the sum of the O₄ optical depth and a low-order polynomial.

Since $\mathrm{\Delta }\stackrel{\mathrm{̃}}{\mathit{\tau }}(\mathit{\lambda },\mathbf{\Omega },N_{\text{a}}(z\left)\right)$ is a non-linear and complicated functional of N_a(z), an analytical solution of this minimization problem does not exist. To simplify the solution of the minimization problem given by Eq. (6), let us consider the variation of DSOT caused by the variation of the aerosol number density. The variation of DSOT in a linear approximation can be represented in the form of the following functional Tailor series (Rozanov and Rozanov, 2010):

$$\begin{array}{ll}\text{(7)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }\stackrel{\mathrm{̃}}{\mathit{\tau }}(\mathit{\lambda },\mathbf{\Omega },N_{\text{a}}(z\left)\right)={\displaystyle }& {\displaystyle }\mathrm{\Delta }\stackrel{\mathrm{̃}}{\mathit{\tau }}(\mathit{\lambda },\mathbf{\Omega },\stackrel{\mathrm{‾}}{N_{\text{a}}}(z\left)\right)+\underset{\mathrm{0}}{\overset{H}{\int }}W(\mathit{\lambda },\mathbf{\Omega },z)\mathit{\delta }{N}_{\text{a}}\left(z\right)\text{d}z+{\mathit{\epsilon }}_{\text{lin}}\left(\mathit{\lambda }\right),\end{array}$$

where

$$W(\mathit{\lambda },\mathbf{\Omega },z)={\left.{\displaystyle \frac{\mathit{\delta }\mathrm{\Delta }\stackrel{\mathrm{̃}}{\mathit{\tau }}(\mathit{\lambda },\mathbf{\Omega },N_{\text{a}}(z\left)\right)}{\mathit{\delta }{N}_{\text{a}}\left(z\right)}}\right|}_{\stackrel{\mathrm{‾}}{N_{\text{a}}}}$$

is a functional derivative of DSOT with respect to the aerosol number density profile N_a(z) around the initial guess $\stackrel{\mathrm{‾}}{N_{\text{a}}}\left(z\right)$. $W(\mathit{\lambda },\mathbf{\Omega },z)$ is usually referred to as the weighting function and ε_lin(λ) is the linearization error.

It follows that the weighting function provides a linear relationship between the variation of DSOT and the variation of an aerosol number density N_a around the initial guess $\stackrel{\mathrm{‾}}{N_{\text{a}}}$. Now inserting Eq. (7) into Eq. (6) and approximating the integral with a finite sum, we have

$$\begin{array}{ll}\text{(8)}& {\displaystyle }& {\displaystyle }\Vert \mathrm{\Delta }\mathit{\tau }(\mathit{\lambda },\mathbf{\Omega })-\mathrm{\Delta }\stackrel{\mathrm{̃}}{\mathit{\tau }}(\mathit{\lambda },\mathbf{\Omega },\stackrel{\mathrm{‾}}{N_{\text{a}}}(z\left)\right){\displaystyle }& {\displaystyle }{-\sum _{k=\mathrm{1}}^L\mathcal{W}(\mathit{\lambda },\mathbf{\Omega },z_k)(x_k-\stackrel{\mathrm{‾}}{x_k})-P(\mathit{\lambda },\mathbf{\Omega })∥}^{\mathrm{2}}⟶min.\end{array}$$

Here, $\mathrm{\Delta }x=x_k-\stackrel{\mathrm{‾}}{x_k}=N_{\text{a}}\left(z_k\right)-\stackrel{\mathrm{‾}}{N_{\text{a}}}\left(z_k\right)$, L is the number of altitude levels, and $\mathcal{W}(\mathit{\lambda },\mathbf{\Omega },z_k)$ is defined to satisfy the trapezoidal integration rule.

The solution of the minimization problem given by Eq. (8) is performed on a discrete wavelength grid. In fact, only five wavelengths over the O₄ absorption band centred, for example, at 477 nm are used.

The minimization problem can be reformulated in the following vector–matrix form:

$$\begin{array}{}\text{(9)}& {∥\mathrm{\Delta }\mathit{y}-\mathbf{K}\mathrm{\Delta }\mathit{x}∥}^{\mathrm{2}}⟶min.\end{array}$$

Here, the vector $\mathrm{\Delta }\mathit{y}=\mathit{y}-\stackrel{\mathrm{‾}}{\mathit{y}}=|\mathrm{\Delta }{\mathit{y}}_{\mathrm{1}},\mathrm{\Delta }{\mathit{y}}_{\mathrm{2}},\mathrm{\dots },\mathrm{\Delta }{\mathit{y}}_N{|}^T$ has the dimension N⋅M and describes the difference between measured and simulated DSOT at N LOS angles and M wavelengths. The mth element of vectors y_n and ${\stackrel{\mathrm{‾}}{\mathit{y}}}_n$ is the differential slant optical depth of O₄ at wavelength m, given by

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }\mathit{\{}{\mathit{y}}_n{\mathit{\}}}_m=\mathrm{\Delta }{\mathit{\tau }}^{-}({\mathit{\lambda }}_m,{\mathbf{\Omega }}_n),m=\mathrm{1},\mathrm{\dots },M,\text{(11)}& {\displaystyle }& {\displaystyle }\mathit{\{}\stackrel{\mathrm{‾}}{{\mathit{y}}_{\mathit{n}}}{\mathit{\}}}_m=\mathrm{\Delta }{\stackrel{\mathrm{̃}}{\mathit{\tau }}}^{-}({\mathit{\lambda }}_m,{\mathbf{\Omega }}_n,\stackrel{\mathrm{‾}}{N_{\text{a}}}\left(z\right)),\end{array}$$

where the superscript “−” denotes that a polynomial is subtracted. The state vector Δx has the dimension L and matrix K has the dimension $(N\cdot M)\times L$; i.e. the matrix consists of N⋅M rows and L columns. The rth row of matrix K is then given by

$$\begin{array}{ll}\text{(12)}& {\displaystyle }& {\displaystyle }\mathit{\{}\mathbf{K}{\mathit{\}}}_r={\displaystyle }& {\displaystyle }\left\{{\mathcal{W}}^{-}({\mathit{\lambda }}_m,{\mathbf{\Omega }}_n,z_{\mathrm{1}}),{\mathcal{W}}^{-}({\mathit{\lambda }}_m,{\mathbf{\Omega }}_n,z_{\mathrm{2}}),\mathrm{\dots },{\mathcal{W}}^{-}({\mathit{\lambda }}_m,{\mathbf{\Omega }}_n,z_L)\right\}\\ {\displaystyle }& {\displaystyle }\text{with}r=(n-\mathrm{1})M+m\end{array}$$

and contains weighting functions for the mth wavelength, nth LOS in L layers.

The minimization problem given by Eq. (9) is solved by employing an iterative Tikhonov regularization technique (Rodgers, 2004):

$$\begin{array}{ll}\text{(13)}& {\displaystyle }{\mathit{x}}_{i+\mathrm{1}}=& {\displaystyle }{\mathit{x}}_{\mathrm{0}}+({\mathbf{K}}_i^T{\mathbf{S}}_y^{-\mathrm{1}}{\mathbf{K}}_i+{\mathbf{S}}_{\mathrm{0}}^{-\mathrm{1}}+\mathit{\gamma }{\mathbf{S}}_t^T{\mathbf{S}}_t{)}^{-\mathrm{1}}{\displaystyle }& {\displaystyle }{\mathbf{K}}_i^T{\mathbf{S}}_y^{-\mathrm{1}}(\mathit{y}-{\mathit{y}}_i+{\mathbf{K}}_i({\mathit{x}}_i-{\mathit{x}}_{\mathrm{0}}\left)\right).\end{array}$$

Here, i is the iteration number, S₀ and S_y are the a priori and measurement covariance matrices, K_i is the weighting function matrix calculated using the estimated aerosol number density profile x_i, S_t is the first-order derivative matrix with Tikhonov parameter γ and ${\mathit{x}}_{\mathrm{0}}=\stackrel{\mathrm{‾}}{\mathit{x}}$ and ${\mathit{y}}_{\mathrm{0}}=\stackrel{\mathrm{‾}}{\mathit{y}}$ are the a priori profile and the representing a priori measurement vector respectively. Equation (13) is solved in the logarithmic space in order to avoid negative values.

There are three criteria to stop the iteration process:

• Convergence in parameter space, i.e. the maximum difference between the components of the state vector at two subsequent iterative steps does not exceed the selected criterion (e.g. 0.0001 km⁻¹);

• The root mean square difference between measured and simulated O₄ DSOT is less than selected (e.g. 0.001);

• The maximum number of iterations is reached.

In addition to the equations above, here we introduce other quantities which are useful for describing the retrieval. The gain matrix

$$\begin{array}{}\text{(14)}& {\mathbf{G}}_i=({\mathbf{K}}_i^T{\mathbf{S}}_y^{-\mathrm{1}}{\mathbf{K}}_i+{\mathbf{S}}_{\mathrm{0}}^{-\mathrm{1}}+\mathit{\gamma }{\mathbf{S}}_t^T{\mathbf{S}}_t{)}^{-\mathrm{1}}{\mathbf{K}}_i^T{\mathbf{S}}_y^{-\mathrm{1}}\end{array}$$

describes the sensitivity of the solution to the measurement. Note that this formulation also includes the Tikhonov term, leading to differences for γ≠0, with the definition given in Rodgers (2004). The averaging kernel A_i=G_iK_i (AK) characterizes the sensitivity of the solution to the true state. The trace of AK quantifies the degrees of freedom of the signal d_s=tr(A) (DOFs). This quantity is commonly understood as the number of individual pieces of information which can be retrieved.

Table 1 SCIATRAN settings for the calculation of differential slant column densities.

3.2 Retrieval of trace gas concentration profiles

Compared to aerosols, the inverse problem for trace gases is easier to solve because, under the assumption of an optically thin atmosphere, the relationship between trace gas concentration and measured differential slant column density is linear. Then, the forward model F(x,b) is equal to a set of measurements y,

$$\begin{array}{}\text{(15)}& \mathit{y}=\mathbf{F}(\mathit{x},\mathit{b})+\mathit{ϵ},\end{array}$$

where ϵ includes the error of measurement and forward model. F(x,b) depends on the retrieval quantity vector x (trace gas concentration profile) and on an additional parameter vector b. The latter includes quantities which have an impact on the measurement and are known with some accuracy (e.g. a priori information).

From the perspective of MAX-DOAS, this relationship is strongly ill posed as the number of retrieval parameters is usually much higher than the number of measurements. Furthermore, as the light paths of two geometrically close elevation angles traverse similar vertical layers near the surface, measurements cannot be considered independent. Consequently, the retrieval parameter vector x is not fully constrained by the input vector y, which introduces the need for additional a priori knowledge which constrains the solution. With the assumption of a Gaussian error distribution, the OE method is often used to prevent unstable solutions when dealing with ill-posed problems.  (Rodgers, 2004).

$$\begin{array}{ll}{\displaystyle }{\mathit{x}}_n& {\displaystyle }={\mathit{x}}_{\mathrm{0}}+({\mathbf{K}}^T{\mathbf{S}}_y^{-\mathrm{1}}\mathbf{K}+g^{-\mathrm{1}}{\mathbf{S}}_{\mathrm{0}}^{-\mathrm{1}}{)}^{-\mathrm{1}}{\mathbf{K}}^T{\mathbf{S}}_y^{-\mathrm{1}}(\mathit{y}-\mathbf{K}{\mathit{x}}_{\mathrm{0}})\\ \text{(16)}& {\displaystyle }& {\displaystyle }\text{or}{\displaystyle }& {\displaystyle }={\mathit{x}}_{\mathrm{0}}+{\mathbf{S}}_{\mathrm{0}}{\mathbf{K}}^T({\mathbf{KS}}_{\mathrm{0}}{\mathbf{K}}^T+g{\mathbf{S}}_y{)}^{-\mathrm{1}}(\mathit{y}-\mathbf{K}{\mathit{x}}_{\mathrm{0}})\end{array}$$

Again, S₀ and S_y are the covariance matrices, y is the measurement and x₀ is the a priori profile. Equation (16) is solved in linear space. Negative values are possible but can be avoided by applying appropriate a priori constraints. Note that we introduced a scaling factor g which gives the possibility of regulating the weighting between a priori and measurement information. Here, the matrix K consists of BAMF values for every layer and measuring geometry (Eq. 4) instead of weighting functions. The covariance matrix of measurements has only diagonal elements, which are the absolute errors of the spectral DOAS fit of this absorber. S₀ has constant variances on its diagonal elements but consists of additional non-diagonal elements based on a Gaussian distribution, which accounts for a possible correlation of different profile layers (Barret et al., 2002).

4 Results and discussion

In this section, we first present and discuss the results of synthetic sensitivity studies. For this purpose, noise-free differential slant columns densities of different aerosol and NO₂ scenarios were simulated with SCIATRAN. Subsequently, these data sets were used as input for BOREAS, and the results are compared with the true profiles. In the next subsection, we discuss the error sources of the retrieval before; in the last section, the results of data measured during the CINDI-2 campaign (2nd Cabauw Intercomparison of Nitrogen Dioxide measuring Instruments, 2016) are shown. These results are validated by comparison with different ancillary measurements (sondes, lidar, long-path DOAS, ceilometer, AERONET, Pandora and in situ).

Table 2 Parameters for aerosol and NO₂ profiles used in the synthetic sensitivity study.

Figure 3 Exponential aerosol profiles for the synthetic sensitivity study. Also shown is the a priori profile for the aerosol retrieval (dashed line).

4.1 Sensitivity study with synthetic data

A synthetic data set of differential slant column densities of NO₂ and O₄ was created using SCIATRAN's spherical mode including multiple scattering. The settings in Table 1 were chosen to describe conditions as they can be found in urban areas like the city of Bremen. Pressure and temperature profiles were taken from a U.S. Standard Atmosphere (NASA, 1976). The vertical grid was set to 25 m steps from the surface up to 6 km, 250 m steps up to 10 km and 1 km steps up to 60 km. The LOS angles were chosen as a representation of a typical MAX-DOAS measurement with more angles for the lower elevations, where the sensitivity is highest. Both the solar zenith angle and relative azimuth angle were kept constant.

Figure 4 Retrieval results with a fixed (a) and a pre-scaled a priori profile (b) with varying Tikhonov parameters γ for SNR = 3000. Small γ values mean less smoothing of the resulting profiles.

The trace gas and aerosol profiles in Table 2 (see also Figs. 3, 6, A1 and A8) were chosen to assess the BOREAS retrieval capability under simple (exponential profiles with scale height SH =1 km) and non-ideal conditions (box profiles, which are not similar to the exponential a priori profiles). The NO₂ scenarios did not include aerosols to avoid the mixture of uncertainties from the trace gas retrieval and inaccuracies within the aerosol retrieval. The box profiles are shown and shortly discussed in Appendix A.

4.1.1 Sensitivity of the aerosol retrieval

Figure 3 shows the true exponential extinction coefficient profiles (see Table 1) with aerosol optical thicknesses of 0.2, 0.6 and 1.0. In addition, an a priori profile is depicted with the same scale height (SH = 1 km) but with AOT = 0.18. In general, an a priori profile should comprise known information of the true atmosphere. Ancillary measurements are frequently used to gather information for a mean profile. Here, we try to evaluate the performance of the algorithm by using an a priori profile close to the true profile (E1) but also for atmospheres with highly varying aerosol loads (E2, E3).

The retrieval uses vertical grid steps of 50 m for synthetic data because smaller step widths introduce retrieval noise, whereas larger steps result in vertical structure being overly smoothed. Since the retrieval tends to overestimate the upper layers due to missing sensitivity (see Fig. 9), a height-dependent a priori variance was used for the preparation of the a priori covariance matrix S_a. The measurement covariance matrix S_e consists of a constant signal-to-noise (SNR) value, which is multiplied with the measurement for each LOS to get absolute covariances. The Tikhonov parameter is changed in the first examples to demonstrate the retrieval response to this parameter. The iterations stop when convergence is reached or after 60 steps. Convergence in parameter space was set to 0.0001 km⁻¹. In addition, convergence of the rms of measured and simulated DSOT differences was set to 0.001.

On the left side of Fig. 4 the retrieval results for different Tikhonov parameters γ are shown. E1, which is close to the a priori profile (dashed grey line) is retrieved well for all γ values. In general, smaller Tikhonov parameters give too much weight to the measurement by reducing smoothing and introducing retrieval noise, which leads to oscillations. Large values smooth the solution in the direction of the a priori profile. When the true atmosphere is far away from the a priori profile, the algorithm does not perform satisfyingly. For E2, strong smoothing (large γ values) leads to more or less straight lines in which insufficient curvature can be found, whereas weak (or no) smoothing means the beginning of oscillations indicated by the reduction of the bottom extinction. No parameter value is able to retrieve the proper bottom extinction. This can be found in the same way for E3, where all results show a clear lack of accuracy. This is not surprising as the a priori profile cannot be understood as the best guess or median of the atmosphere any more.

Table 3 True and retrieved AOT for SNR = 3000 and γ=1.0 with and without a priori pre-scaling, including relative differences to the true profile in brackets.

Figure 5 Variation of the relative difference in AOT and bottom extinction (BOT) compared to the true value as a function of SNR and Tikhonov parameter, colour coded for scenario E3 with a priori pre-scaling.

Several authors have already highlighted retrieval problems from MAX-DOAS measurements under highly variable (in time and space) atmospheric conditions (Frieß et al., 2018; Vlemmix et al., 2011, 2015; Wagner et al., 2011). Nevertheless, studies that focus on improvements of the regularization (Eriksson, 2000; Hartl and Wenig, 2013) or a priori information (Clémer et al., 2010; Frieß et al., 2006; Hendrick et al., 2014) are rare with respect to MAX-DOAS retrievals. Here, we use a priori pre-scaling to cope with strong deviations of the true atmosphere from the a priori estimate. For this purpose, the BOREAS aerosol retrieval is started on a reduced vertical and spectral grid from a zero profile (all a priori extinction values were set to be zero) as a priori information. In Eq. (15), the a priori profile x₀ is set to 0 but the remaining terms stay unchanged in order to retrieve a reliable AOT value without limiting the solution by a certain profile shape. Since the retrieval works directly on O₄ absorption, aerosol optical thickness is in general the most reliable retrieval quantity. The AOT from this pre-calculation is used to scale the a priori profile for the main run by calculating a scaling factor from the initial a priori AOT and the new value, which is applied to all initial extinction coefficients before the main inversion. Note that the only a priori information used both within the pre-run and within the main run is the Tikhonov parameter. This method performs well when the a priori shape is similar to the true atmospheric profile shape (right side of Fig. 4). However, no improvement is found when the shape of the a priori profile deviates strongly from the true atmospheric condition but no deterioration is found either (see box profiles in Appendix A). On the right-hand side of Fig. 4, clear improvements can be seen in both the overall shape and bottom extinction. Nevertheless, the profiles tend to deviate more strongly at the surface, especially when the aerosol load is higher, which was also found by Frieß et al. (2006). The relative differences in retrieved AOT are high for the unscaled retrieval but can be significantly reduced with a priori pre-scaling (see Table 3). Further improvements could be made by increasing the allowed number of iterations together with stronger convergence criteria (not shown), but the average retrieval time would increase as well.

A comparison of all three profiles and their optimal γ ranges with and without pre-scaling indicates that each combination of a priori and true profile has a specific optimal Tikhonov parameter which might differ strongly between each scenario. In addition, the SNR is another quantity which has a large impact on the retrieval. In general, SNR values might differ strongly throughout a day, because dynamic ranges of exposure and integration times are frequently used to operate the CCD in the ideal saturation range when the number of photons might differ strongly (e.g. due to different viewing geometries). In Fig. 5 we varied SNR (as a representation of the measurement error) and the Tikhonov parameter γ from 500 to 5000 and from 0.25 to 30 respectively to find the best parameter range for the optimization of AOT and the bottom extinction (BOT). The figure shows results for the strong aerosol scenario E3 (E1, E2 and the box scenarios are shown in the appendix). The relative differences to the true value are colour coded with bright spots, indicating better retrievals.

For scenario E3, high γ (or low SNR) values apply too much smoothing (not enough weight on the measurement), which ends in strong AOT ∕ BOT underestimations of the retrieval results. Furthermore, there is an ideal range of SNR values for γ≤2 where the AOT relative difference is ≤10 %. In this specific SNR ∕ γ range, the bottom values do not show the best results. On the other hand, BOT values are optimized for $\mathrm{3}\le \mathit{\gamma }\le \mathrm{10}$ for SNR =2500 and SNR =3000 where the AOT does not have the lowest difference. As a consequence, due to unequal vertical sensitivity, the choice of a priori constraints and parameterization, in some cases either AOT or the shape of the profile can be optimized. Stable solutions for the optimization of both profile characteristics could not be identified in the range of parameterization, regularization and convergence. As a result, a reasonable compromise for optimum values with respect to both AOT and BOT would be in the range SNR = 2000–3000 with low Tikhonov values $\mathrm{2}\le \mathit{\gamma }\le \mathrm{5}$. Note that the SNR range between 3500 and 4000 has larger relative deviations than the surrounding settings which could not be found for other scenarios. This is assumed to be either a unique solution, where the range of settings leads to poor retrieval results, or a numerical instability.

In contrast to this pragmatic approach, the so-called L-curve approach (Engl and Grever, 1994; Hansen, 1992) was found to be insufficient as it already assumes good regularization; i.e. the data residual and the parameter norm (∥x−x_a∥) can be described well by a specific regularization parameter (Hartl and Wenig, 2013), which was not found for the application of MAX-DOAS profiling.

4.1.2 Sensitivity of the trace gas retrieval

The optimal-estimation-based trace gas profile retrieval faces similar problems to the regularization difficulties discussed for the aerosol retrieval in the previous section. In many studies, an a priori variance of a fixed percentage of the a priori profile is used in combination with the differential slant column density errors from the DOAS fit to constrain the solution (Frieß et al., 2018). Here, we show that this is an insufficient approach, since the fitting errors might be too small, giving too much weight to the measurement when the NO₂ concentration is high and the a priori profile is not a good estimate of the true atmosphere.

Again, a priori pre-scaling is used to achieve a better estimate of the true atmospheric conditions. In contrast to the aerosol retrieval, which needed a pre-run for a better first guess, the trace gas a priori pre-scaling is achieved using the 30^∘ elevation angle measurement. In the geometric approximation (see Eq. 3), the differential AMF is 1, which leads to ΔS₃₀=V. This vertical column density V is used as a pre-scaling value for the a priori profile.

The different NO₂ exponential profiles, based on the parameters in Table 2, are shown in Fig. 6 together with the a priori profile (box profiles can be found in Appendix A). Note that we used an a priori profile close to E2. These profiles are a good representation of the daily variability in an urban region, where high NO₂ concentrations are found during phases of commuter traffic (rush hours) in the morning and late afternoon and low NO₂ loads are found in between.

Figure 6 Exponential NO₂ profiles for the synthetic sensitivity study. Also shown is the a priori profile for the aerosol retrieval (dashed line).

Figure 7 depicts the profiling results without a priori pre-scaling on the left side and with pre-scaling on the right side. Different regularization ratios are achieved by varying parameter g in Eq. (16), which varies the measurement error by $\sqrt{g}$. The retrieval is done on the same vertical grid as in the previous section with similar variances for the a priori covariance matrix. In contrast to a constant SNR, we use variances in the range of typical DOAS fit errors of NO₂ for the measurement covariance matrix, which can be understood as the standard approach of the MAX-DOAS profiling community. For small elevation angles and high NO₂ concentrations, the relative ΔS error is often lower than 1 % and sometimes even lower than 0.3 %, indicating a high reliability in the measurement but neglecting errors due to pointing inaccuracies of the telescope or horizontal inhomogeneities which affect the profiling result but not the measurement.

The results again show the importance of pre-scaling. Whereas the E2 scenario is well retrieved for all g ranges without pre-scaling, the other scenarios tend to oscillate or underestimate the profile. Note that the standard regularization ratio (g=1) does not lead to satisfying results for E3. Applying pre-scaling leads to more accurate profiles, especially for E3.

Figure 7 Retrieval results with a fixed (a) and a pre-scaled a priori profile (b) with g factors. Small g values mean less measurement weighting for the resulting profiles.

Figure 8 Retrieval errors for E1–E3 of the synthetic aerosol (a) and NO₂ study (b). Aerosol results are shown for SNR = 3000 and Tikhonov parameter γ=2. Trace gas errors are depicted without additional regularization (g=1). Pre-scaling was used.

The specific optimal g range for each scenario depends not only on the a priori profile and the true atmosphere but also on the variance values, their height dependencies and the definition of off-diagonal elements. However, the quantification of this optimal range is a difficult task for real data due to the lack of knowledge about the true atmosphere. Manual variations of g as well as iterative approaches are conceivable. Possible criteria for finding the best g values are, for example, the rms between measurement and retrieval quantity or ancillary measurements. Here, we used fixed values which were chosen to decrease oscillations while maintaining the best surface concentration.

4.2 Discussion of error sources and information content

The total error of both aerosol and trace gas retrieval can be separated into three different error sources (Rodgers, 2004).

$$\begin{array}{}\text{(17)}& {\displaystyle }{\mathbf{S}}_{\text{tot}}={\mathbf{S}}_{\text{sm}}+{\mathbf{S}}_{\text{fw}}+{\mathbf{S}}_{\text{ns}}\end{array}$$

The first error term describes the smoothing error ${\mathbf{S}}_{\text{sm}}=(\mathbf{A}-\mathbf{I}){\mathbf{S}}_{\text{e}}(\mathbf{A}-\mathbf{I}{)}^T$ which can be understood as error for the averaging kernel smoothed estimate of the true state of the atmosphere. A and I are the averaging kernel and identity matrix respectively. The second term of S_tot is the forward model error S_fw, which is hard to quantify in real retrievals, because it needs the true state of the atmosphere. For near-linear problems, this error source can be neglected, when the forward model parameters are well estimated. The last error is called the retrieval noise S_ns=GS_eG^T (G: gain matrix; see e.g. Eq. 14) and denotes the uncertainty of the retrieved profile due to measurement errors.

Figure 8 shows smoothing, retrieval and total error for the above-mentioned scenarios of the aerosol retrieval (left) and the trace gas retrieval (right).

The total errors for both aerosol and trace gas retrieval are dominated by the smoothing error with negligibly small retrieval noise. This shows that the measurement itself can be considered a small error source in contrast to the generally limited information content of MAX-DOAS measurements. Note that the error components of the NO₂ profile show small oscillations, which also indicates that larger constraints are needed. For the aerosol retrieval, the largest errors can be found at an altitude at which the sensitivity is low, with an additional increase near the surface. The errors are nearly zero for higher altitudes because of the height-dependent variances, which were chosen in order to reduce possible retrieval instabilities in altitude regimes where the sensitivity is lowest. In general, NO₂ also shows larger errors for higher altitudes but the surface error is dominant here, indicating a reduced sensitivity or poor parameterization. The errors for z≥3 km are close to zero because the smallest a priori covariance is used as a boundary condition for the highest altitudes, where the sensitivity is low.

Figure 9 Averaging kernels of the aerosol retrievals for scenarios E1–E3 (first three figures on the left side). Also shown are the area and FWHM of the AK on the right side. FWHM values are depicted on their nominal heights (coloured) and on the height of their individual peaks (black).

For the quantification of the vertical sensitivity of the retrieval, Fig. 9 shows the averaging kernels AK (see Eq. 14 and text) for the aerosol scenarios E1–E3 calculated for a Tikhonov parameter γ=0.25 and SNR = 3000. We used a small γ to demonstrate the general vertical sensitivity rather than the smoothed one. Larger γ values would lead to peaks at lower altitudes with slightly reduced peak sizes. Within an ideal retrieval, the AK matrix should be equal to the identity matrix (Rodgers, 2004). Then, delta peaks in every layer would indicate a perfect sensitivity of the retrieval to the true state of the atmosphere. For MAX-DOAS retrievals, AK are much broader, with peaks close to the surface. Here, the sensitivity is highest and the kernels for a certain layer show correlations with neighbouring layers due to the finite width of each kernel. The comparison of AK for E1–E3 reveals that, besides the surface kernel, all peaks are reduced in size and width for higher aerosol loads. This demonstrates the smaller sensitivity due to both an insufficient a priori estimate and the additional aerosol load, which reduces the general number of photons scattered into the light path (for specific geometries, e.g. RAA =90^∘). The area of the AK should ideally be 1 over a large range of altitudes. Here, we find that only within the first kilometre is the overall sensitivity high. Note that, especially for E2 and E3, the surface area value is much smaller than 1 due to the negative area of the lowermost kernels for higher altitudes. Negative kernel values indicate that an additional aerosol load in this altitude would decrease the retrieved solution close to the surface. The full width at half maximum (FWHM) of averaging kernels are a measure of the vertical retrieval resolution with the lowest values close to the surface. In the right-hand-side plot of Fig. 9 we show FWHM at the nominal altitudes (coloured lines) and the altitudes of their individual maxima (black lines). The calculation was done by using the individual AK instead of the FWHM of fitted Gaussian distributions. The limited vertical extent of the black curves indicate that a theoretical vertical resolution with values between 500 m and 1 km is available for all kernel levels, but the specific sensitivity for all layers lies within the lowest kilometre.

The trace of the averaging kernel matrix (degrees of freedom, DOFs) is often discussed as if it is the number of independent pieces of information for a retrieval. Hartl and Wenig (2013) already pointed out that this number is not necessarily a strict limitation of the retrieval and might differ strongly depending on the regularization.

In Fig. 10, we depicted the rms between measured and retrieved O₄ ΔS and the DOFs for scenario E3. Giving the retrieval more freedom will increase the degrees of freedom and decrease the rms, which might be understood as more information is used for a better measurement assimilation. However, in Sect. 4.3.1, we already demonstrated that profiles with high DOFs and low rms might oscillate. DOFs only state the maximum number of independent pieces of information for the current regularization parameters without considering if the retrieved profiles are optimal. Therefore, we do not use this quantity for the description of retrieval results.

Figure 10 Rms between measured and simulated O₄ ΔS and the matching degrees of freedom (DOFs) for the specific retrievals for different SNR and γ values.

In Fig. 11, we show the averaging kernel, area and FWHM of the NO₂ retrieval for scenarios E1–E3. Similar conclusions can be drawn as for the aerosol retrieval for the shape and altitude variation of AK. The area is closer to 1 for all altitudes, which is due to the missing smoothing of the Tikhonov term. The FWHM again shows a lower vertical sensitivity for higher layers. On the other hand, the FWHM curves at their specific peak heights (black) show a larger vertical spread in comparison to the aerosol curves, which indicates a larger sensitivity for higher altitudes in the NO₂ retrieval. The DOFs for the trace gas retrieval lie in the range between 2.3 and 3.2 but decrease significantly when applying values of g≠1.

Figure 11 Averaging kernels of the NO₂ retrievals for scenarios E1–E3 (first three figures on the left side). Also shown are the area and FWHM of the AK on the right side. FWHM values are depicted at their nominal heights (coloured) and at the heights of their individual peaks (black).

In addition to the above mentioned errors and resolution properties, further error sources and problems should be noted. For real measurements, the true meteorological conditions are mostly unknown. Wang et al. (2017) pointed out that wrong O₄ profiles might lead to strong deviations between retrieved and true aerosol profiles. Furthermore, we did not implement the instruments field of view (FOV) in the retrieval because it will increase run time. However, larger differences can be expected for concentrations and extinctions in the lowermost layers, when parts of the FOV might gather photons reflected on the Earth's surface. Additionally, we assume a linearity between measurement and concentration, which might be inaccurate for high aerosol loads and NO₂ concentrations when the telescope points to the horizon.

4.3 Retrieval of profiles during the CINDI-2 campaign

The 2nd Cabauw Intercomparison of Nitrogen Dioxide Measuring Instruments campaign (CINDI-2) took place in Cabauw (the Netherlands) from 25 August to 7 October 2016 and was funded by the European Space Agency (ESA) and the participating research groups. The campaign goals were the characterization of differences between NO₂ measurements by various instruments and approaches, and the evaluation of the resulting data sets for the validation of the Copernicus satellite, Sentinel 5 Precursor (S5P). Over 40 different instruments operated by 30 groups from all over the world provided an outstanding ensemble of data sets for this task. The campaign was a successor of CINDI, which was held in Cabauw from 16 to 24 June 2009 (Pinardi et al., 2013; Piters et al., 2012; Zieger et al., 2011).

The measurement site Cabauw is located in a rural region dominated by agriculture but is surrounded by four of the largest cities in the Netherlands (Rotterdam, Amsterdam, Den Haag and Utrecht). Thus, depending on the wind direction, long-range transport from highly industrialized areas is likely, which results in high-pollution events. Here, 3 of the 5 investigated days show a more or less steady wind from south-easterly directions, with a change in wind direction in the evening of 15 September (see Fig. 12). On 23 September, the wind came from the west, whereas 1 day later the wind came mostly from southerly directions.

Figure 12 Wind speed and direction for the investigated days (CESAR, 2018) resampled to 5 min values from 06:00 to 17:00 UTC.

In this study, the instrumental set-up for the validation of aerosol and trace gas profiles consists of in situ and remote-sensing instruments. Near-surface concentration and extinction values are provided by a ceilometer, NO₂ lidar, long-path DOAS (LP DOAS) and in situ samplers (ICAD, CAPS, NAQMN). Integrated values are provided by Pandora and AERONET direct Sun instruments. Besides the MAX-DOAS measurements, the only routinely retrieved profiling information comes from ceilometer and NO₂ lidar data. In addition, three NO₂ sondes were launched on one of the days investigated here (15 September 2016).

The profile validation is realized on 2 cloud-free days (13–14 September 2016) and 3 days with broken clouds (15 September and 23–24 September 2016). The operators of MAX-DOAS instruments were asked to perform elevation scans at the beginning of each hour and at 11:15 and 11:45 UTC every day. The ancillary measurements introduced in the following subsections were resampled or averaged in this time interval to prevent discrepancies due to time lags between two observations. Each scan took around 8 min and consisted of the following elevation angles: 1, 2, 3, 4, 5, 6, 8, 15 and 30^∘. The telescope pointed towards the westerly direction with a viewing azimuth angle of 107^∘ from the south.

In this section, BOREAS retrievals were performed on data from the IUPB MAX-DOAS instrument on a 100 m step width vertical grid from the surface up to 4 km. The aerosol retrieval uses the a priori pre-scaling option with an exponential a priori profile described by a 1 km scale height and an aerosol optical thickness of 0.18. The asymmetry factor and single-scattering albedo were chosen to be 0.92 and 0.68 respectively. The SNR was set to 2500 with a Tikhonov parameter of 2. The a priori variance decreased with altitude from 1.5 at the surface to 0.01 at 4 km. This definition was chosen to improve the profiling results for lower altitudes, where the main aerosol load/concentration can be expected, and in order to suppress instabilities at the upper grid boundary, where the sensitivity is lowest. The trace gas retrieval again uses a priori pre-scaling with an exponentially decreasing profile with a scale height of 1 km and a vertical column density of 9×10¹⁵ molecules cm⁻². The a priori covariance matrix uses the same variances as the aerosol retrieval but includes Gaussian distributed side diagonal elements to account for correlations between individual layers (Barret et al., 2002). The measurement covariance matrix for NO₂ consists of the total differential slant column error on the diagonal elements only. Vertical profiles of pressure and temperature were created by taking the mean of 16 different sonde measurements taken during the years 2013–2015 in De Bilt (the Netherlands).

Figure 13 (a) Time series of AERONET (blue) and BOREAS AOT (red). Small triangles show the original AERONET measurement, small dots with connecting lines depict the resampled data. Grey areas indicate clouds. (b) Scatter plot of both data sets including parameters of the orthogonal regression and Pearson's correlation coefficient.

Figure 14 (a) Time series of BOREAS (red), ceilometer (green) and NAQMN in situ (blue) and near-surface aerosol parameters. Small triangles show the original hourly NAQMN measurements; small dots with connecting lines depict the resampled data. Green squares show ceilometer near-surface extinction values evaluated by averaging the 10–50 m coefficients. Grey areas indicate clouds. (b) Scatter plot including parameters of the orthogonal regression and Pearson's correlation coefficients.

4.3.1 Validation of aerosol retrievals

Instrumentation for the aerosol validation

Aerosol profiles were validated using three different instruments. The retrieved AOT was compared with values from an AERONET station (AErosol RObotic NETwork, Dubovik et al., 2000; Holben et al., 1998). The level 2 AOT at 440 nm was scaled with the Angström exponent from the ratio 440 nm/675 nm to calculate the AOT at 477 nm (see Fig. 13). Due to measurement intervals varying between 4 and 30 min, we decided to resample the AERONET signal. The errors from AERONET instruments are usually in the range of 0.01 (VIS, IR) to 0.02 (UV) (Sayer et al., 2013). Here, we used a constant error of 0.01.

The bottom extinction coefficient of the retrieved aerosol profiles was compared with in situ PM₁₀ concentrations from the National Air Quality Monitoring Network (NAQMN) operated by the Dutch National Institute for Public Health and the Environment (RIVM). The NAQMN measurement site Wielsekade is located at a distance of around 900 m from the MAX-DOAS station. The site is listed as a regional background station and the PM₁₀ measurement principle is beta attenuation (Thermo Fisher Scientific FH62I-R). The instrument provides data on 1 min intervals but only the hourly values are validated by RIVM. These hourly values were resampled on the MAX-DOAS times and the minute values were used to calculate the errors as standard deviation for the data points shown (see Fig. 14).

Furthermore, we used AERONET-scaled ceilometer near-surface extinction as a further validation data set. The ceilometer (CHM15k Nimbus) was operated by the Royal Netherlands Meteorological Institute (KNMI) and sampled backscattering signals every hour at a wavelength of 1064 nm. The integrated backscattering signal was divided by the 1020 nm AERONET AOD to get a conversion factor which was applied to the backscattering signal for the conversion into extinction coefficients. This new ceilometer profile was again scaled with AERONET AOT at 477 nm. The error is expressed as the standard deviation calculated from the temporal and vertical averaging.

Comparison of aerosol retrieval parameters

Figure 13 shows the time series and the scatter plot of AERONET- and BOREAS-derived AOT. The BOREAS data were filtered for profiles which reached the iteration limit during the retrieval. Grey areas indicate clouds within the specified time period.

In the scatter plot, BOREAS error bars were calculated by integrating the total error (Eq. 17) vertically. On the first 2 cloud-free days, the temporal variability shows a similar pattern with a slight offset between the instruments, most likely introduced by an underestimation of the AOT by BOREAS. Due to the limited sensitivity of MAX-DOAS profiling for higher altitudes (see Fig. 9), elevated aerosols do not contribute to BOREAS AOT but can be measured by direct Sun instruments like AERONET Sun photometers.

In the evening, the results of both instruments vary more. Different reasons for this finding are possible. First, a developing planetary boundary layer (PBL) might exceed the vertical extension where BOREAS has sufficient sensitivity. Second, when pointing towards the Sun, saturation of the CCD becomes a problem, leading to low integration times, which decrease SNR and fitting quality. In addition, RTM calculations might introduce uncertainties at high SZA, when the aerosol load is high and the light path is strongly increased. Furthermore, the aerosol phase function used leads to uncertainties, as the forward scattering peak might be underestimated by the Henyey–Greenstein parameterization. The other days show a higher variability due to an increase in cloudy scenes. Note that, especially in the morning and around noon, clouds may influence the instruments in various ways because of their different azimuthal viewing directions. The BOREAS profile retrieval cannot be considered reliable when one or several elevation angles point at clouds during the measurements, as intensity and light path vary in unpredictable ways. This might also lead to large iteration numbers or to no convergence at all. On 15 September, the temporal patterns of the two instruments differ strongly. Although the extinction values around noon agree well, AERONET's AOT decrease before and the increase afterwards cannot be found this strongly by BOREAS. Especially in the evening, both curves deviate, indicating highly variable atmospheric conditions, which might also be introduced by a change in wind direction (see Fig. 12). The aerosol load on the last 2 days is much lower because of rainfall in the time between 16 and 22 September. In the morning of 23 September, several profiles were discarded due to exceeded iteration limits.

The correlation coefficient of 0.83 is high and shows the general good agreement on 3 of 5 days. The regression line indicates the above-mentioned small underestimation of the AOT by BOREAS.

Figure 14 shows the comparison of near-surface aerosol values. Extinction coefficients from BOREAS and the ceilometer are depicted together with PM₁₀ concentrations of the NAQMN in situ instrument. Note that the y axis for PM₁₀ concentrations was chosen to match the results of BOREAS on 14 September. The NAQMN in situ sampler and BOREAS show a good temporal agreement within the first 3 days. Only in the morning hours of the 13 and 15 September, do data from the two instruments show a larger spread. This deviation is most likely due to the insufficient vertical resolution of MAX-DOAS profile retrievals when the PBL has not yet evolved and near-surface aerosol loads are dominant (see Figs. 9 and 15 for further details). On 15 September, early morning clouds might also have an impact on the deviation. Again, a larger variability for BOREAS results can be found for the last 2 days, when the retrieval was influenced by broken clouds. The in situ instruments show a smoother daily variation here. The ceilometer bottom extinction shows a good agreement with BOREAS and NAQMN on the first day and an underestimation on the second day. On 15 September, bottom values are influenced by thick clouds which might interfere with the backscattering signal of altitude ranges below the clouds (see Fig. 16). The day 23 September shows more variability than the first 2 days but this can be found for the other instruments as well.

The correlation of BOREAS bottom extinctions with data from both validation instruments is high, indicating good agreement on cloud-free days. The correlation with ceilometer near-surface values is high as well but the regression line shows a general underestimation of ceilometer near-surface values in comparison to BOREAS.

Figure 15 shows the aerosol layer height found by the ceilometer for all days. Especially during the first 3 days, an upward-extending PBL can be identified from noon to the late afternoon. Only on 15 September, does a more or less stable residual PBL seems to exist in the morning hours. The other days show high variability in the morning with layer heights from 200 to 600 m with more or less individual high signals, which might be produced by smaller clouds at higher altitudes: 24 September shows a stable layer around 200 m and a diffuse developing PBL in the afternoon.

Figure 15 Ceilometer aerosol layer height within the planetary boundary layer.

The underlying extinction coefficient profiles for the above-discussed bottom extinction and AOT values can be seen in Fig. 16. In the top row, temporal and vertically averaged ceilometer profiles are depicted. In the middle row, these data were smoothed to the MAX-DOAS vertical resolution by the application of BOREAS averaging kernels (${\mathit{x}}_{\text{new}}={\mathit{x}}_{\text{apri}}+\mathbf{A}({\mathit{x}}_{\text{ceilo}}-{\mathit{x}}_{\text{apri}}$); see Rodgers and Connor, 2003). The first day shows a good agreement between BOREAS and ceilometer near-surface extinction with a small offset between the curves. In the morning, both instruments find the aerosol load mainly located in the lowermost layers. Beginning around noon, the PBL starts developing with a maximum height of 2 km in the afternoon found by the ceilometer. BOREAS cannot resolve this increasing PBL as well as the averaged ceilometer but the AK smoothed data indicate that the reason might be limited sensitivity to the top boundary of the PBL.

On 15 September, the averaged ceilometer data show high and thick clouds in the morning and a rising PBL in the afternoon. In the AK-smoothed data, these clouds cannot be identified any more. BOREAS introduces elevated aerosol layers which can be understood as a retrieval artefact due to these cloudy scenes. In the afternoon, BOREAS finds an upward expansion in the PBL similar to the ceilometer, with the exception of the last profile, which was influenced by the telescope pointing towards the Sun. BOREAS extinction values are smaller in the PBL, which is a consequence of the already explained underestimation of BOREAS AOT to AERONETs AOT, which was used for the backscatter signal scaling.

Figure 16 Comparison of aerosol extinction coefficient profiles from BOREAS and ceilometer for 13 September 2016 (a) and 15 September 2016 (b).

4.3.2 Validation of nitrogen dioxide profiles

Instrumentation for the NO₂ validation

Several instruments were used for the validation of BOREAS nitrogen dioxide profiling results. The VCD results were validated with the help of a Pandora instrument (no. 128) operated by LuftBlick (Herman et al., 2009) and a NO₂ lidar from RIVM. In addition, three NO₂ sondes launched on 15 September by KNMI were also used.

The Pandora instrument was used with its direct Sun capability to retrieve total nitrogen dioxide columns on hourly time steps (see Fig. 17). Since the measurement time did not match to BOREAS scans, Pandora values were resampled. The stratospheric column was subtracted based on the approach of Knepp et al. (2015) by using OMI stratospheric columns during the Cabauw overpass (L2 OVP, NASA, 2018). The data were quality filtered and the errors were provided by LuftBlick and are based on Herman et al. (2009).

The NO₂ lidar (Brinksma et al., 2008) did not measure as frequently as Pandora and MAX-DOAS instruments but performed several measurements on 15 September. One lidar scan took approximately 30 min with the temporal midpoints of the scans on a non-regular temporal grid which necessitated resampling. The scan was done on different elevation angles which enabled the retrieval of vertical profiles with altitude points up to 2.5 km. For the VCD calculation, these profiles were vertically integrated. Errors were provided by the RIVM team and Gaussian error propagation was applied for the calculation of VCD errors. In addition to the remote-sensing instrument, NO₂ sondes were launched on 15 September around 05:15, 08:04 and 10:25 UTC. The sondes provide NO₂ VMR as well as pressure, temperature and altitude data for the conversion to concentrations, which were vertically integrated up to 4 km.

Near-surface concentration validation was done with the help of in situ (NAQMN, ICAD, CAPS) and remote-sensing instruments (lidar and LP DOAS). The in situ samplers NAQMN, CAPS and ICAD were operated by RIVM, the Royal Belgian Institute for Space Aeronomy (BIRA-IASB) and the Institute of Environmental Physics in Heidelberg (IUPH) respectively. In addition to the aerosol bottom concentration, NAQMN also provided NO₂ concentrations (Teledyne API 200E; see Sect. 4.3.1 for details on data handling). IUPH and BIRA-IASB created a common NO₂ in situ data set from IUPHs ICAD (Iterative CAvity DOAS, Pöhler et al., 2017) and BIRAs CAPS (Environment SA, AS32M-CAPS, Cavity Attenuated Phase Shift, Kebabian et al., 2008). This joint data set filled gaps in the individual measurement series and fixed a scaling issue for CAPS. BIRA-IASB operated two different CAPS instruments, which were installed at the 27 and 200 m levels of the Cabauw measurement tower. The joint ICAD/CAPS data set as well as the 200 m CAPS data were averaged over 15 min time steps. For both instruments (ICAD/CAPS) an error of 1 ppb was assumed.

Figure 17 (a) Time series of Pandora, lidar and BOREAS NO₂ VCD. Blue diamonds show the original Pandora measurements and green squares show the lidar data points. Small dots with connecting lines depict the resampled data sets of Pandora and lidar. BOREAS is shown as red circles. Grey areas indicate clouds. The plot includes orange triangles as a representation of the integrated NO₂ sondes measurements with the ascent and descent separated into different triangles with the edge at the top or bottom respectively. (b) Scatter plot of the data sets including parameters of the orthogonal regression and Pearson's correlation coefficients.

Figure 18 (a) Time series of in situ and remote-sensing NO₂ near-surface concentrations. Triangles show in situ instruments for NAQMN (bright blue, triangles pointing upwards) and for ICAD/CAPS (dark blue, triangles pointing downwards). In addition, on 13 September 2016, CAPS at the 200 m level is depicted as green triangles with the edge up and the mean values of both CAPS with the edge to the side (yellow). The NO₂ lidar is plotted as green squares and the LP DOAS as different-sized circles in shades of magenta (lowest altitude as smallest circle). BOREAS is shown as large red circles. (b) Scatter plot of the data sets including parameters of the orthogonal regression and Pearson's correlation coefficients.

In addition to the in situ instruments, a long-path DOAS (LP DOAS) provided by IUPH (Pöhler et al., 2010) was used for the validation. Four reflectors attached at different altitude levels of the measurement tower (12.7, 47, 107 and 207 m) provided unique and well-defined light paths which enabled the calculation of concentration on several altitudes. The instruments' elevation measurements led to profiles every 30 min. Errors were calculated by the IUPH team.

Nitrogen dioxide profiles were provided only by the NO₂ lidar and sonde launches. The LP DOAS did not cover enough altitude levels for a comparison. The conversion from trace gas VMR to concentrations for all data sets was done with the meteorological data from the CESAR observatory (CESAR, 2018).

Comparison of NO₂ retrieval parameters

Figure 17 depicts the comparison of VCD for different instruments. The agreement with Pandora is good on all days with larger deviations in the morning and on 15 September in general with the exception around noon. Deviations in the morning might be (similarly to the AERONET comparison) due to spatial differences in the NO₂ distribution, when Pandora points to the east and MAX-DOAS to the west. This argument is supported by the lidar measurement on 15 September that agreed well with BOREAS and shared the same azimuthal viewing direction. In the afternoon of that day, BOREAS VCD increased faster than Pandora and lidar, with a strong overestimation for the last profile at 16:00 when the Sun was low. Here, we can exclude horizontal gradients as a reason for the differences, because of the same viewing direction of lidar and MAX-DOAS. On the same day, NO₂ sondes proved that the BOREAS VCD was correct only during the ascent (triangle with edge at the top) of the second launch. The matching descent (triangle with edge at the bottom) for the first flight was approximately 38 km in the north-westerly direction between Amsterdam and Rotterdam and agrees more with Pandora columns, indicating differences due to azimuthal viewing directions. The strong differences in VCD for ascent and descent of the first and second sonde launches show high temporal and spatial variability of the NO₂ concentration, which is supported by all remote-sensing instruments.

The correlation with Pandora is slightly better than with lidar but the regression parameter shows that Pandora VCDs are higher than those from BOREAS, whereas the lidar regression is closer to 1.

In Fig. 18, near-surface NO₂ concentrations are depicted. Both in situ data sets (NAQMN and ICAD/CAPS) agree well with each other and show larger differences only in the morning of 14 September. The LP DOAS profiles at 12.7 and 47 m show similar concentrations to the in situ instruments. Differences in the morning of the first 2 days between the lower and the higher LP DOAS values indicate a strong vertical inhomogeneity, which was also found for the aerosols (see Sect. 4.3.1).

BOREAS surface concentrations agree very well with all data sets except during the morning hours. On 13 September, in addition to the ICAD/CAPS data at the 27 m level (blue), the 200 m points are shown (green). The anti-correlated behaviour in the morning hours confirms the strong inhomogeneity and proves that nitrogen dioxide was mainly concentrated close to the surface at that time. The mean values of both CAPS instruments (yellow line) show a better agreement with BOREAS than the 27 and 200 m data.

It seems that BOREAS cannot fully resolve this thin near-surface layer and retrieves rather smooth profiles instead of sharp concentration peaks (see also the discussion about vertical sensitivity in Sect. 4.2). This is again in agreement with the LP DOAS results on 24 September, where the 200 m level LP DOAS finds similar concentrations to BOREAS. Note that the NO₂ lidar has a finer resolution within the lowest 100 m, which might be the reason it depicts similar concentration numbers to the other instruments in the morning hours of 15 September. As a rough estimation of the thickness of the individual thin layers in the morning, we divided BOREAS VCD by the ICAD/CAPS concentration and found that the layer height lies around 200 m (06:00–08:00), 300 m (09:00) and 400 m (10:00) on 13 September. The mean AK FWHM for the surface layer on that day is 214 m, indicating that the lowermost layer can be resolved only from around 09:00 UTC where the curves approach each other. A good agreement is found at 10:00 UTC when the concentration is focused in a layer twice as thick as BOREAS surface resolution.

The correlations are high for all instruments with the highest value for NAQMN and the NO₂ lidar. The in situ instruments show similar correlations, but a slight underestimation of BOREAS can be found with slopes between 1.22 and 1.44.

Figure 19 depicts NO₂ profiles for BOREAS, lidar and sondes on 15 September. Similarly to the aerosol profile comparison, averaging kernels were applied to lidar and sondes measurements. The vertical profiles of the unsmoothed lidar show the previously discussed high vertical inhomogeneity. We can clearly identify a high concentration close to the surface and another elevated NO₂ layer with altitudes varying between 100 and 500 m. The unsmoothed sondes agree well with the lidar measurements only at the ascent of the second launch. After the application of averaging kernels, the two distinct layers with high concentrations are smoothed to one layer with NO₂ concentrated at the surface, again indicating a lower vertical sensitivity of MAX-DOAS profiling. In addition, smoothed and unsmoothed lidar find larger concentrations for altitudes over 500 m before 08:00 UTC, which was not found by BOREAS and the sonde profiles.

Figure 19 Comparison of NO₂ concentration profiles from BOREAS, lidar and sonde measurements for 15 September 2016.

5 Conclusions

In this study, we introduced BOREAS, a new profiling algorithm for MAX-DOAS measurements. BOREAS retrieves aerosol extinction coefficient profiles around O₄ absorption bands from the oxygen dimers differential slant column densities. In contrast to existing inversion algorithms, which are based on perturbation theory, BOREAS directly minimizes the difference between modelled and measured optical depth by varying aerosol extinction coefficient profiles in an iterative Tikhonov regularization scheme. In addition, trace gas concentration profiles are retrieved via optimal estimation using the aerosol properties determined in the first step. For both retrieval parts, we introduced a priori pre-scaling as a simple method for improving results under variable atmospheric conditions.

The retrieval performance was demonstrated on the example of synthetic data calculated with SCIATRAN and with real measurements taken during the CINDI-2 campaign in Cabauw 2016. The synthetic scenarios were chosen to work as a stress test for BOREAS. One fixed a priori profile was used for the retrieval of different large aerosol loads and NO₂ concentrations. The results are good when the scenario is relatively close to the a priori profile but fails for strong loads and concentrations. In a comprehensive variation of parameters we showed that, in non-ideal situations, either the shape/near-surface value or the integrated quantity can be optimized. The information content of a MAX-DOAS scan with degrees of freedom between ≈1 and ≈3.5 is not enough for a perfect retrieval. Furthermore, the frequently utilized assumption of using the trace gas errors from DOAS fits only as measurement variance was found to be insufficient since the regularization ratio might be wrongly weighted. We suggest a manual change in measurement errors or an iterative approach for finding the best ratio.

In the second part of this study, BOREAS retrieval results from measurements of the IUPB MAX-DOAS instrument during the CINDI-2 campaign were validated with ancillary measurements. The agreement for all parameters was good with correlations equal or higher than 0.75. Systematic offsets or deviations depending on geometry, atmospheric condition and investigated air mass indicate limitations due to the specific measurement characteristics, which were also found in earlier MAX-DOAS profiling studies. The correlation with AERONET AOT is 0.83 with an underestimation by BOREAS, leading to a regression slope of 1.30. The bottom value correlation is 0.75 for the ceilometer near-surface extinction and 0.81 for NAQMN in situ measurements. The total column of NO₂ shows a good agreement with Pandora and lidar with correlations of 0.86 and 0.79 respectively. In contrast to the underestimation of BOREAS integrated aerosol, Pandora nitrogen dioxide columns are overestimated, especially in the morning hours. The correlations for near-surface NO₂ concentrations are high with values larger than 0.84. Stronger deviations between BOREAS and in situ instruments were found in the morning. This was explained by the limited vertical resolution of BOREAS, while in the morning the highest concentrations are found close to the surface.

As a conclusion, BOREAS profiling capabilities are strong, which was proven by high correlations with all validating instruments. Discrepancies were found due to different azimuthal viewing directions of the instruments and the limited vertical resolution of MAX-DOAS profiling when aerosol load or concentration is close to the surface in layers thinner than the BOREAS vertical resolution. A comparison of BOREAS results to the performance of other MAX-DOAS profiling algorithms will be discussed in two upcoming papers (Frieß et al., 2018; Tirpitz et al., 2019).

Appendix A: Additional plots and retrieval of box profiles

A1 Additional aerosol scenarios and plots

Figure A1 Box aerosol profiles for the synthetic sensitivity study. Also shown is the a priori profile for the aerosol retrieval (dashed line).

The retrieval results for the box scenarios in Fig. A1 are shown in Fig. A2 for different Tikhonov parameters and SNR = 3000. A priori pre-scaling leads to a clear improvement only for scenario B3 since the bottom extinction is reached. For B2 and B3, the resulting profile shapes change slightly, which indicates the above-mentioned lack of improvements through pre-scaling, when the a priori profile shape differs strongly from the true atmosphere. Nevertheless, since exponential profiles are improved a lot and box profiles show no degradation, we highly recommend the usage of better and more flexible a priori profiles.

Figure A2 Retrieval results with a fixed (a) and a pre-scaled a priori profile (b) with varying Tikhonov parameters γ for SNR = 3000. Small γ values mean less smoothing of the resulting profiles.

Figure A3 Variation of SNR and Tikhonov parameter with the relative differences in AOT and bottom extinction compared to the true value, colour coded for scenario B1 with a priori pre-scaling.

Figure A4 Same as Fig. A3 but for scenario B2.

Figure A5 Same as Fig. A3 but for scenario B3.

Figures A3–A5 show the SNR and Tikhonov variation for the three aerosol box profiles. One can see that for B2 the smallest γ values lead to the best bottom value, whereas for scenario B1 these values would result in an overestimation. Also the optimal SNR ∕ γ range for the AOT differs strongly for all box profiles, which again shows that regularization should be considered a flexible choice of parameters, depending on the specific atmospheric conditions.

Figures A6 and A7 show the SNR and Tikhonov variation for the exponential profiles E1 and E2. Note that for E1 all parameter combinations lead to good profiling results because the a priori profile is close to the true atmosphere. In comparison to E3, E2 shows a better adaptation of AOT with a coincidently similar performance for the bottom extinction.

Figure A6 Same as Fig. 5 but for scenario E1.

Figure A7 Same as Fig. 5 but for scenario E2.

Additional NO₂ scenarios and plots

Figure A8 shows the box NO₂ scenarios, which are retrieved for different g values in Fig. A9. The standard regularization ratio for g=1 is again found to be oscillating without reaching the proper bottom concentrations for E2 and E3. Higher g factors show an improvement for both the fixed a priori and pre-scaled initial guess.

Figure A8 Box NO₂ profiles for the synthetic sensitivity study. Also shown is the a priori profile for the aerosol retrieval (dashed line).

Figure A9 Retrieval results with a fixed (a) and a pre-scaled a priori profile (b) with g factors. Small g values mean less measurement weighting for the resulting profiles.