Retrieval of tropospheric aerosol, NO 2 and HCHO vertical proﬁles from MAX-DOAS observations over Thessaloniki, Greece: Intercomparison and validation of two inversion algorithms

. In this study we focus on the retrieval of aerosol and trace gas vertical proﬁles from Multi-Axis Differential Optical Absorption Spectroscopy (MAX-DOAS) observations for the ﬁrst time over Thessaloniki, Greece. We use two independent inversion algorithms for the proﬁle retrievals: The Mexican MAX-DOAS Fit (MMF) and the Mainz Proﬁle Algorithm (MAPA). The former is based on the Optimal Estimation Method (OEM), while the latter follows a parameterization approach. We evaluate the performance of MMF and MAPA and we validate their retrieved products with ancillary data measured by other co-5 located reference instruments. The trace gas differential slant column densities (dSCDs), simulated by the forward models, are in good agreement, except for HCHO, where larger scatter is observed due to the increased spectral noise of the measurements in the UV. We ﬁnd an excellent agreement between the tropospheric column densities of NO 2 retrieved by MMF and MAPA (Slope = 1 . 009 , Pearson’s correlation coefﬁcient R = 0 . 982 ) and a good correlation for the case of HCHO ( R = 0 . 927 ). For aerosols, we ﬁnd better agreement for the aerosol optical depths (AODs) in the visible (i.e., at 477 nm), compared to the UV (at 10 360 nm) and we show that the agreement strongly depends on the O 4 scaling factor that is used in the analysis. The agreement for NO 2 and HCHO near-surface concentrations is similar to the comparison of the integrated columns with slightly decreased correlation coefﬁcients. The seasonal mean vertical proﬁles that are retrieved by MMF and MAPA are intercompared and the seasonal variation of all species along with possible sources is discussed. The AODs retrieved by the MAX-DOAS are validated by comparing them with AOD values measured by a CIMEL sun-photometer and a Brewer spectrophotometer. Four

∼ 1 o . Simultaneous azimuth and elevation angle calibration is regularly performed by sighting the sun, so no horizon scans are necessary for the elevation angle calibration. A routine MAX-DOAS measurement cycle starts by orienting the optics at a certain azimuth viewing direction followed 110 by the measurement of scattered radiation spectra at the elevation angles: 90 (zenith), 30, 15, 12, 10, 8, 6, 5, 4, 3, 2 and 1 o in this order. For this study, the system was configured to measure at four consecutive azimuth angles of 142, 185, 220 and 255 o , illustrated in Figure 2 with arrows of different colors. Based on the intensity of the measured spectra during an elevation scan at 142 o azimuth, the viewing direction of 1 o elevation angle was found to be partly blocked by obstacles, such as trees and buildings in the campus. Thus, α = 1 o in this particular direction was excluded from the profiling analysis. In order to 115 achieve high SNR values and to avoid saturated spectra, the number of scans of each individual measurement and the exposure time of the CCD are automatically adjusted by the operating software according to the received intensity by the detector. The integration time at each elevation angle is ∼ 60 sec and a full measurement sequence for all azimuth directions lasts about one hour.  Vandaele et al. (1998), I0-corrected (SCD in 10 17 molecules cm -2 ) NO2 (220 K) Vandaele et al. (1998) (2000) H2O (296 K) HITEMP, Rothman et al. (2010) Ring Ring spectra calculated by QDOAS according to Chance and Spurr (1997) Polynomial degree 5 5 5 Intensity offset Constant Order 1 Constant Wavelength Calibration Based on a high-resolution solar reference spectrum (Chance and Kurucz, 2010) 2.4 MMF 165 The Mexican MAX-DOAS Fit (MMF) v2020_04 (Friedrich et al., 2019) is an OEM-based profiling algorithm that relies on online RTM simulations using VLIDORT version 2.7 (Spurr, 2006) as forward model. The input parameters for each atmospheric layer are calculated from temperature and pressure profiles, the trace gas concentration in each layer and the aerosol properties. The aerosol properties, which are the same for all layers, are the single scattering albedo (SSA) and the asymmetry parameter (using the Henyey-Greenstein phase function, Henyey and Greenstein (1941), to calculate the phase 170 function moments). Furthermore, the wavelength of the retrieval and the surface albedo need to be specified as additional input parameters. The retrieval algorithm comprises an aerosol extinction profile retrieval and a trace gas profile retrieval. The former constrains the aerosol extinction profile in the forward model of the trace gas retrieval. The inversion uses constrained damped least-square fitting with an optimal estimation regularization. In the used version, both the a priori and the covariance matrix are constructed. More details about the a priori settings and the input parameters can be found in Sect. 2.6. The retrieval algorithm  (Beirle 195 et al., 2019), which can describe the majority of all potential measurement sites. MMF, on the other hand, relies on online RTM simulations and so the aerosol and surface parameters can be adjusted to the most suitable values. In this study, the aerosol optical properties that are used as input for the simulations of MMF are based on 15 years climatological data measured by a co-located CIMEL sun-photometer. Figure 4 shows the frequency distribution of the Ångström exponent, AOD, asymmetry factor and SSA in Thessaloniki, while their values that are used as input to each inversion algorithm are listed in Table 2.

200
Discrepancies between MMF and MAPA due to small differences in these selected parameters are expected to be minor. MMF requires a priori profile and covariance matrix information for the profile retrievals. The "a priori" term represents knowledge of the true state before the measurement is performed. However, the true shape of the trace gas vertical profiles at Thessaloniki is generally not known, while the true state of the aerosol profiles is known only in certain cases during the period of study. Thus, the retrieval is based on constructed exponentially decreasing a priori profiles with scale height of 1 km, which 205 are considered a reasonable estimate of the true profiles. Since no covariance matrix information is available, the covariance matrix is also constructed from the a priori profile. The AOD as well as the trace gas VCDs in Thessaloniki vary substantially throughout the year. In order to take into account the annual variability, we use the square of 50% of the a priori on the diagonal elements of the covariance matrix for aerosols and 100% for NO 2 and HCHO. The loose constraint of the latter is due to the higher variability of the trace gas vertical columns over the course of the year. Both for aerosols and trace gases, the off-axis 210 elements of the covariance matrix were constructed by assuming a correlation length of 200 m, as described in Clémer et al.
(2010). Additionally, the progress of the convergence is faster when using an a priori VCD or AOD below the true value. Thus, the a priori AODs were set to 0.25 and 0.15 for the aerosol retrievals at 360 and 477 nm, respectively. For the trace gas retrievals we have used a priori VCDs of 4 × 10 15 and 6 × 10 15 molecules cm -2 for NO 2 and HCHO, respectively, based on data derived from the MAX-DOAS by applying the geometrical approximation (Hönninger et al., 2004) to the dSCDs measured at 30 o and and 0.2 -1.8 for the profile shape parameter (Beirle et al., 2019).
The temperature and pressure vertical profiles that are used as input in this study are identical for both MMF and MAPA.
We have used climatological profiles for Thessaloniki produced by MPIC that are based on ∼ 16 years re-analysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF). The temperature and pressure profiles are interpolated 220 to the day and time of each elevation sequence. We have also tried to use temperature and pressure profiles measured by radiosondes, launched on a daily basis at Thessaloniki Airport (∼ 13 km away from the measurement site), as input, but since no major effect is observed on the retrieved products, these results are not presented. Both algorithms are configured to export the retrieved vertical profiles to the same output grid ranging from the ground up to 4 km with 200 m vertical resolution. Aerosol Robotic Network (AERONET) (https://aeronet.gsfc.nasa.gov/). CIMEL is an automated, well-calibrated scanning filter radiometer specifically developed for the retrieval of the AOD at 7 wavelengths (i.e., 340, 380, 440, 500, 670, 870 and 1020 nm) by using direct-sun observations. The technical specifications of the instrument are given in Holben et al. (1998).

230
The instrument is calibrated regularly following the procedures and the guidelines of AERONET. The AERONET database provides three distinct levels for data quality. Level 1.0 is defined as pre-screened data (i.e., no quality assurance criteria are applied). The Version 3 (Sinyuk et al., 2020) of Level 1.5 represents near-real-time automatic cloud screened data, while Level 2.0 applies additional pre-and post-field calibrations. In this paper, we use the AERONET Level 1.5 data, since the Level 2.0 data for the period of study is not yet published. In order to compare with the AOD retrieved by the MAX-DOAS, the AODs 235 at 360 and 477 nm have been calculated using the Ångström exponent between the standard spectral bands of the instrument.

Brewer spectrophotometer
The Brewer spectrophotometer with serial number 086 (B086) is a double monochromator that performs spectrally resolved measurements of the direct and global solar irradiance at Thessaloniki since 1993 (Bais et al., 1996;Fountoulakis et al., 2016).
The wavelength range of B086 is 290 -365 nm and its spectral resolution is 0.55 nm at full width at half maximum (FWHM).

240
The wavelength calibration is performed by scanning the emission lines of spectral discharge lamps, while maintenance of the absolute calibration is achieved by regularly scanning the spectral irradiance of a calibrated 1000-W quartz-halogen tungsten lamp (Garane et al., 2006).
Although the Brewer's initial purpose was the retrieval of total ozone columns, research activities have shown that the spectral AOD can be calculated from direct irradiance measurements by following two main approaches: The first is based on 245 the absolute calibration of the direct-sun spectra measured by the Brewer (Kazadzis et al., 2005), while the second uses the Langley extrapolation method (relative calibration) (Gröbner and Meleti, 2004). In both cases the spectral AOD is calculated as the residual optical depth after subtracting from the total atmospheric optical depth the optical depths due to molecular scattering and the O 3 and SO 2 absorption (Kazadzis et al., 2007). Since 1997, the direct solar irradiance spectra measured by the Brewer are calibrated (Bais, 1997), so in this study we use the former approach (i.e., absolute calibration) for the retrieval 250 of the spectral AOD. In order to compare with the AOD retrieved by the MAX-DOAS, the AOD at 477 nm is calculated using climatological monthly mean values of the extinction Ångström exponent derived from measurements of the CIMEL sun-photometer in Thessaloniki. Details on the procedure of Brewer's direct solar irradiance spectra absolute calibration, as well as the spectral AOD retrieval methodology can be found in Bais (1997); Kazadzis et al. (2005Kazadzis et al. ( , 2007; Fountoulakis et al. (2019).

Lidar
Thessaloniki is a member station of the European Lidar Aerosol Network (EARLINET, https://www.earlinet.org) since 2000, providing regular aerosol profile measurements, following the EARLINET's schedule (Monday morning, Monday and Thursday evening), during extreme events and at satellite overpasses (e.g., AEOLUS, OMI).
THEssaloniki LIdar SYStem (THELISYS) is a multi-wavelength Raman/depolarization lidar system, which has been grad-260 ually upgraded regarding its operational wavelengths and the detection configuration. All the quality standards, established within EARLINET, are followed in order to assure the high quality of the THELISYS products, which are publicly available in the EARLINET database (https://www.earlinet.org/index.php?id=125). A detailed description of THELISYS technical specifications and algorithm can be found in Voudouri et al. (2020).
The final products derived from the raw lidar data processing are: the aerosol backscatter coefficient at 355, 532 and 1064 nm, 265 the aerosol extinction coefficient at 355 and 532 nm and the linear particle/volume depolarization ratio at 532 nm. During the day, the data acquisition is limited to the signals that arise from the elastic scattering of the laser beam by the air molecules and the atmospheric aerosol. The Klett-Fernald algorithm in backward integration mode is applied (Klett, 1981)  Another source of uncertainty during the lidar signals processing is the system's overlap function, which determines the altitude, above which a profile contains trustworthy values. In our analysis, the correction is not available for the daytime 275 retrievals. Thus, an overlap function from the previous nighttime measurement or a mean overlap profile is applied. The starting height is set to the full overlap height (approximately 0.6 km), assuming height-independent backscatter below 0.6 km, equal to the backscatter measured at this height, to account for both the incomplete overlap within the lidar profile and atmospheric variability in the lowermost tropospheric part. This overlap effect generally introduces uncertainties in the calculation of the columnar products (e.g., AOD unaffected by local traffic emissions and therefore can be considered more representative of the average NO 2 concentrations in the local boundary layer.

Results and discussion
In this section, we present results of the trace gas and aerosol quantities retrieved by the two inversion algorithms. We intercompare the integrated columns (i.e., VCDs and AODs for trace gases and aerosols, respectively), the dSCDs simulated 295 by the RTMs, the retrieved vertical profiles and the surface concentrations between MMF and MAPA. Since MAPA is based on a parameterization approach, no information about averaging kernels is provided; hence, results on averaging kernels are presented only for MMF.
The MAX-DOAS system operates at a site where the northern viewing directions are blocked by buildings of the campus and the city, so the system is configured to perform sequences of elevation scans at azimuth directions in the southern sector, 300 as illustrated in Figure 2. As a result, scattered radiation spectra may be measured during the day at azimuths close to the solar azimuth angle. In such cases, RTM simulations might face difficulties in calculating properly the dAMF due to increased aerosol forward scattering, leading usually to underestimation of the true dAMF. Therefore, the elevation sequences measured at azimuth angles relative to the sun of less than 5 o are excluded from the analysis. In addition, the elevation sequences, for which the retrieved AOD from the MAX-DOAS inversion algorithms is greater than 1.5 are filtered-out, since such high 305 aerosol loads are unrealistic for Thessaloniki ( Figure 4). Negative columns can occur in the trace gas retrievals of MAPA within the Monte Carlo ensemble and they are intentionally not removed, but this is not possible for MMF retrievals since, in its current version, MMF operates in logarithmic state vector space. For NO 2 , no valid negative columns are retrieved, but for HCHO, MAPA reports negative columns for ∼ 8.5% of the valid data. In order to compare meaningful results between the two algorithms, the negative columns are removed from the initial dataset. It should be noted that in the following sections an Orthogonal Distance Regression (ODR or bivariate least-squares) has been used instead of an Ordinary Linear Regression (OLR or standard least-squares) for the comparison of the retrieved 325 products derived by MMF and MAPA, in order to treat equally the two algorithms since none of them depends on the other.
The discrepancies in the regression slopes and intercepts arising in the OLR when comparing independent variables, and the appropriateness of ODR, are discussed in Cantrell (2008). The ODR results are also sensitive to the assumed errors of the two variables. The uncertainty contained in the MAX-DOAS measurements may be difficult to assess, but, since both MMF and MAPA retrievals are based on the same input data, the associated errors are assumed the same and equal to the mean error 330 provided by MMF and MAPA for each data point.

Integrated Columns
In the past, the trace gas VCDs measured by our MAX-DOAS systems have been derived by dividing the measured dSCDs, only at two elevation angles, 30 o and 15 o , or the mean of the two, with appropriate dAMFs. The dAMFs have been calculated either following the geometrical approximation approach or by deploying RTM simulations taking into account the viewing 335 geometry, the aerosol optical properties and the instrument's viewing direction relative to the sun (Drosoglou et al., 2017).
However, in both cases, the actual trace gas profile has not been taken into consideration introducing, possibly, an additional uncertainty to the measured VCD. This is the first time during the Phaethon's operation that the whole elevation sequence is being used in order to derive the tropospheric VCDs more accurately.
In Figure 5 the time series of the integrated columns of all retrieved species (i.e., AODs for aerosols and VCDs for trace  for O 4 VIS). Additionally, the simulated dSCDs are in good agreement with the measured dSCDs, which is a good indicator for successful profile retrievals. In the case of O 4 (UV), even though the slope and correlation coefficient are similar to O 4 (VIS), a larger scatter is evident, while for HCHO larger deviations from unity in the slopes and correlation coefficients are observed, especially at higher elevation angles. This can probably be explained by the increased noise in the UV spectra compared to the VIS range and also due to the fact that at higher elevation angles the measured differential optical densities are very low, 375 reaching the spectrometer's detection limit. For aerosols in both spectral ranges, discrepancies between the simulated dSCDs of MMF and MAPA may also arise due to the variable O 4 scaling factor that is included in MAPA retrievals.

Surface concentrations
The surface concentration is defined as the trace gas amount at ground level. However, in MAPA's configuration the algorithm allows for the retrieval of lifted trace gas layers for a shape parameter greater than 1, which leads to a zero value for the surface 380 concentration. For these cases, the comparison of the lowermost concentrations with in situ measurements or surface concentrations retrieved by an OEM-based algorithm will be biased. Thus, in the following sections, the term "surface concentration" will refer to the concentration reported for the lowermost layer of the MAX-DOAS profile for both MMF and MAPA, rather than the concentration at the ground. Figure

Averaging Kernels
The averaging kernels (AVKs) of a profile retrieval describe the sensitivity of the retrieved state to the true atmospheric state for each altitude layer. The degrees of freedom (DoF) are mathematically derived as the trace (or sum of the diagonal elements) of the AVK matrix and quantify the number of independent pieces of information gained from the measurements compared to the a priori knowledge (Rodgers, 2000). Both the AVKs and the DoF can be used to characterize the quality of the retrieved profile.

395
Since only OEM-based inversion algorithms are capable of providing AVKs, the results shown here are derived only by MMF. Figure 8 shows a typical example of the calculated AVKs for each of the retrieved species. The DoF of this example retrieval are shown for each species. The median DoF retrieved by MMF are: 3.13 ± 0.32 for NO 2 , 2.22 ± 0.34 for HCHO, 2.73 ± 0.28 for aerosols in the VIS and 2.02 ± 0.32 for aerosols in the UV. The averaging kernels verify that MAX-DOAS measurements are typically less sensitive for altitudes greater than ∼ 2 km, as a result of the viewing geometry, and thus, altitudes greater than 400 3 km are not presented here. That means that the MAX-DOAS measurements under these viewing geometries and with the a priori profiles and covariance matrices used in this study (Sect. 2.6) are adequate for retrieving the extinction and concentration profiles only up to the lowermost ∼ 1.5 -2 km of the atmosphere with highest sensitivity closer to the ground. Also, since the photon path increases with wavelength, the MAX-DOAS technique shows higher sensitivity for the species retrieved in the VIS range than in the UV.

Validation
In this section we present the validation results of the products retrieved by the MAX-DOAS profile analysis against ancillary data measured by other reference co-located instruments. Vertical profiles of the aerosol extinction measured by a co-located raised by MMF and MAPA should be treated as valid or not. Figure 10 shows the comparison between the common AOD data derived by CIMEL, Brewer and the MAX-DOAS at 360 and 477 nm. Each column of the figure corresponds to a different flagging scheme as described in Table 4. Figure 11 presents graphically the statistics of the linear regressions (i.e., slope, offset, number of points and Pearson's correlation coefficient) between the reference instruments and the MAX-DOAS. The panels a -d correspond to different flagging schemes, as Figure 10.

430
The comparison results of the MAX-DOAS against the CIMEL are slightly different than those with the Brewer. This can probably be explained by the fact that only few collocated measurements are available and in different periods for the two reference instruments (Figure 9). In the case of the AOD at 477 nm most of the outliers are filtered-out when the flagging scheme #4 is applied and the best agreement is observed between the reference instruments and the MAX-DOAS, both for  The effect of the O 4 scaling factor is also investigated by comparing the integrated columns of MMF and MAPA and also by comparing the AODs derived by MAPA for different values of the scaling factor with AODs measured by the CIMEL and the Brewer. The effect of the O 4 scaling factor has a stronger impact on aerosols than on trace gases (where the effect is minor). The fixed value of 0.8 for the scaling factor, which is supported by many studies, does not seem to be suitable for the measurements 555 at Thessaloniki.

Appendix A: Effect of the O 4 scaling factor
The O 4 scaling factor (SF) was introduced by Wagner et al. (2009) in order to remove the systematic discrepancies appearing between measured and simulated O 4 dSCDs. Uncertainties of the O 4 cross sections and/or its temperature and pressure dependence, aerosol optical properties and RTM errors have been suggested as possible causes for these discrepancies. Several 560 studies have confirmed the idea of the O 4 scaling factor and have shown that applying a SF (commonly using a value between 0.75 and 0.9) is indeed necessary (Wagner et al., 2009;Clémer et al., 2010;Irie et al., 2011;Vlemmix et al., 2015b;Wang et al., 2016;Frieß et al., 2016;Wagner et al., 2019Wagner et al., , 2021. However, other studies have not supported this requirement (Spinei et al., 2015;Ortega et al., 2016;Seyler et al., 2017;Wang et al., 2017a, b). Although the need for an O 4 scaling factor for retrieving aerosol information from MAX-DOAS measurements has been extensively discussed (Wagner 565 et al., 2019), its physical mechanism is not understood and still remains an unresolved issue.
Since MAPA provides the option of scaling the modeled O 4 dSCDs, we have investigated the effect of the SF on the comparisons between the products of the two profiling algorithms and between AOD derived by MAPA and the reference instruments. We selected three fixed values (i.e., 1, 0.9 and 0.8, referred hereafter as SF1.0, SF0.9 and SF0.8) and a variable SF (SFvar). Figure A1 shows the effect of the O 4 SF on the comparison of trace gas VCDs and AOD derived by MAPA and 570 MMF (same as Figure 5 without accounting for the different azimuth directions). Like in Figure  The AERONET Version 3 aerosol retrieval algorithm, associated uncertainties and comparisons to Version 2, Atmospheric Measurement Techniques, 13, 3375-3411, https://doi.org/10.5194/amt-13-3375-2020, 2020.