Impact of 3D Cloud Structures on the Atmospheric Trace Gas Products from UV-VIS Sounders-Part III: bias estimate using synthetic and observational data

Three-dimensional (3D) cloud structures may impact atmospheric trace gas products from ultraviolet-visible (UVVIS) sounders. We used synthetic and observational data to identify and quantify possible cloud-related bias in NO2 tropospheric vertical column densities (TVCD). The synthetic data were based on high-resolution large eddy simulations which were input to a 3D radiative transfer model. The simulated visible spectra for low-earth orbiting and geostationary geometries were analysed with standard retrieval methods and cloud correction schemes that are employed in operational NO2 satellite 5 products. For the observational data the NO2 products from the TROPOspheric Monitoring Instrument (TROPOMI) were used while the Visible Infrared Imaging Radiometer Suite (VIIRS) provided high spatial resolution cloud and radiance data. Cloud shadow fraction, cloud top height, cloud optical depth, solar zenith and viewing angles, were identified as the metrics being the most important in identifying 3D cloud impacts on NO2 TVCD retrievals. For a solar zenith angle less than about 40◦ the synthetic data show that the NO2 TVCD bias is typically below 10%. For larger solar zenith angles both synthetic and 10 observational data often show NO2 TVCD bias on the order of tens of %. Specifically, for clearly identified cloud shadow bands in the observational data, the NO2 TVCD appears low-biased when the cloud shadow fraction > 0.0 compared to when the cloud shadow fraction is zero. For solar zenith angles between 50-60◦, about 16% of TROPOMI pixels with high quality value NO2 TVCD retrievals, were found to be impacted by cloud effects larger than 20%.

Partnership (S-NPP) satellite. The S-NPP payload includes the Visible Infrared Imaging Radiometer Suite (VIIRS) instrument which may be used as an imager for TROPOMI. We utilized data from both sensors.
The paper first discusses the synthetic and observational data sets, section 2. This is followed by a description of the various metrics used to identify cloud impacts, section 3. In section 4 the results are presented and discussed before the paper ends by some concluding remarks and an outlook, section 5.

2 Data
The study area covers approximately Germany, the Netherlands and parts of other surrounding countries. For this region we have generated synthetic satellite data using high resolution LES data as input for 3D radiative transfer simulations.

Synthetic satellite data
The synthetic satellite data comprise results from three sources: 1) high-resolution LES cloud data; 2) 3D radiative transfer 10 modelling of satellite radiances with the LES cloud data as input; and 3) NO 2 retrieval using the synthetic satellite radiances.
Within the HD(CP)2 project LES were made with the ICOsahedral Non-hydrostatic atmosphere model (ICON;Dipankar et al., 2015;Zängl et al., 2015) which simulated realistic liquid and ice clouds. The results have a spatial resolution of approximately 1.2×1.2 km 2 for a region including Germany, the Netherlands and parts of other surrounding countries and was validated against ground and satellite-based observational data by Heinze et al. (2017). Several weeks of simulations are avail-15 able including all kinds of weather situations in Europe. We, however, utilize only one from 29 July 2014, 12:00 UTC, due to the computational burden of the radiative transfer simulations. The simulated scene includes all cloud types that are typical for Europe, such as shallow cumulus, cirrus, stratus, and also convective clouds. While focus is on Europe, the results/methods are expected to be general and thus applicable elsewhere.
The LES results were input to the 3D MYSTIC Monte Carlo radiative transfer model (Mayer, 2009;Emde et al., 2011) run 20 within the libRadtran package (Mayer and Kylling, 2005;Emde et al., 2016) to generate synthetic observation spectra in the visible spectral range from 400-500 nm and in the O 2 -A band region from 755-775 nm (for further details see: Emde et al., 2021;Yu et al., 2021). The spatial resolution of the simulated sensor was set to approximately 7×7 km, corresponding to 98×104 pixels for the full LES domain. Note, that each simulated sensor pixel includes 36 cloud pixels, hence the simulations include sub-pixel cloud inhomogeneity. The synthetic spectra were input to a NO 2 retrieval algorithm which included two steps: first 25 a differential optical absorption spectroscopy (DOAS) fit was performed using the QDOAS retrieval algorithm (Blond et al., 2007;De Smedt et al., 2008) to get the NO 2 slant column densities; second the slant column densities were converted to vertical column densities using layer air mass factors based on the VLIDORT 1D radiative transport model (Spurr, 2006). The fitting window was between 425 and 495 nm, similar to the one used by Richter et al. (2011). The air mass factor was calculated at the middle of the fitting window, that is 460 nm. Cloud corrections were made using both O 2 -O 2 and O 2 -A band based 30 methods. An example of retrieved NO 2 TVCD for a synthetic case is shown in Fig. 1a. Note that the "true" NO 2 is constant over the scene with column density of 1.6×10 16 molec/cm 2 corresponding to an European tropospheric polluted NO 2 profile from Levelt et al. (2009). Thus any differences between the retrieved and "true" NO 2 TVCDs are due to the presence of clouds, While we discuss our synthetic results in connection with observational results from satellites in low-earth-orbit (LEO), simulations were also made for a geostationary orbit (GEO). Azimuth and zenith viewing and solar angles were chosen to resemble geometries for the study region when viewed by the TROPOspheric Monitoring Instrument (TROPOMI, Veefkind 5 et al., 2012) and the future Ultra-violet Visible Near-infrared (UVN, https://sentinel.esa.int/web/sentinel/missions/sentinel-4) instrument to be in geostationary orbit. In total 15 and 36 combinations of viewing and solar angles were simulated for the GEO-case and LEO-case, respectively (Emde et al., 2021). The surface was assumed to be snow free and with constant albedo to simplify the interpretation of the results. In the visible (400-500 nm) simulations were made with albedos of 0,0.05 and 0.2, while in the O 2 -A band region an additional albedo of 0.5 was included to account for the potentially larger albedo in this part The TROPOMI data includes the latitude and longitude of the corners of each TROPOMI pixel. This information was used to identify VIIRS pixels within each TROPOMI pixel.
Data within the study region roughly covering Germany and parts of neighbouring countries were collected for the periods 28 June 2018 -15 October 2018 and 1 March 2019 -30 June 2019 when the sun was sufficiently high on the horizon to avoid problems with NO 2 retrievals for low sun (solar zenith angle < 60 • ). An example of data from near simultaneous TROPOMI 5 and VIIRS overpasses of Germany is shown in Fig. 2.

Methods
Several metrics were calculated to quantify cloud features and their possible connection with NO 2 -biases due to 3D cloud effects.

Cloud geometric and radiance fractions
Cloud fractions may be defined in several ways. We calculate the geometric cloud fraction, CF g , the radiometric cloud fraction, CF r , and the weighted radiometric cloud fraction, CF w .
The geometric cloud fraction for a TROPOMI pixel is defined as where the sum is over all N VIIRS pixels within the TROPOMI pixel and CM i is the VIIRS cloud mask where CM i = 1 for pixels identified as cloudy and CM i = 0 otherwise.
The radiometric cloud fraction is the fraction of measured radiance reflected from clouds in a pixel (see also Grzegorski et al., 2006) Here R is the observed reflectance, R s is the reflectance for a cloudless sky and R c the reflectance for an opaque cloud. For the O 2 -O 2 cloud correction, CF r is calculated based on the reflectance at 460nm which is in the middle of the DOAS fitting window for NO 2 . For the FRESCO algorithm, CF r is determined by the reflectance in the 758-759nm window band. Further details are described by Yu et al. (2021). We also define an average radiometric cloud fraction CF V IIRS r using the average of CF r calculated for each VIIRS M3 band pixel, centred at 0.488 µ, within a TROPOMI pixel: where the sum is over all N VIIRS pixels within the TROPOMI pixel.
Finally, the weighted radiometric cloud fraction is defined as (Yu et al., 2021)

Cloud mask 20
No cloud mask was available for the LES based synthetic data. In principle this may be calculated from the cloud optical thickness assuming some treshold. However, we chose to calculate a cloud mask similar to how it is done for several operational satellite products. We thus assume that pixels with reflectance larger than some threshold, R > R t , are cloudy. The cloud mask was calculated from high spatial resolution radiance simulations at 0.55 µm of the synthetic cases described by Emde et al. (2021). A threshold of R t = 0.25, similar to Heinze et al. (2017), was adopted. It is noted that Barker et al. (2017) used 25 R t = 0.15 for the GOES-13 0.65 µm band. Yang and Di Girolamo (2008) has shown that there may be overlap between clear and cloudy pixels and hence some mis-classification is unavoidable. Our slightly higher value potentially includes more cloud contaminated pixels, but also reflects that we are mainly working over land where surface albedo is higher than over ocean.

Cloud shadow fraction and cloud shadow index
The VIIRS cloud shadow mask algorithm is geometry-based and described by Hutchison et al. (2009). They compared the spectral procedure used in the MODIS MOD35 product with the geometry-based approach and found the latter to be far superior. The cloud shadow fraction, CSF , for a TROPOMI pixel is defined as where the sum is over all N VIIRS pixels within the TROPOMI pixel and CSM i is the VIIRS cloud shadow mask where CSM i = 1 for pixels identified as cloud shadow and CSM i = 0 otherwise.
For the LES we define the cloud shadow index (CSI) as a surrogate for a cloud shadow product as follows: Here E 0 = 1 is the direct transmittance at the surface for an atmosphere with no molecular absorption nor clouds or aerosol.

10
The direct transmittance at the surface for a solar zenith angle of θ and solar azimuth angle φ, E θ,φ , was calculcated with the MYSTIC model. It also does not include molecular absorption nor aerosol, however, it includes clouds, which are fully accounted for in 3D geometry and with the 3D cloud fields in higher spatial resolution than the TROPOMI pixel size. The CSI is zero for pixels with no cloud shadow and 1 if a pixel is fully in the cloud shadow. The CSI clearly and unambigously identify cloud shadow regions as it is based solely on sun, cloud and surface geometry. An example of the CSI is given in 15 Fig. 1c.

H-metric
For a TROPOMI pixel the H-metric is the standard deviation of VIIRS band reflectances divided by the mean of the VIIRS reflectances within the TROPOMI pixel. It is an estimate of the variation of radiance within the TROPOMI pixel and has been used earlier by for example Massie et al. (2017) as metric for the impact of 3D clouds on CO 2 retrievals. In the presence of 20 sub-pixel clouds the H-metric is expected to increase as the horizontal cloud inhomogeneity increase. In the absence of clouds the H-metric is expected to be small. However, for cloud free pixels with large variations in surface albedo the H-metric may still be large, and conversely, for a completely cloudy pixel the H-metric may be small. Thus, it may not be used without care to unambiguously identify cloud inhomogeneity. From the VIIRS data we calculate the H-metric for the M3, M4 and M5 bands centred at 0.488, 0.555 and 0.672 µm respectively.

Absorbing aerosol index
Clouds have an effect on the UV absorbing aerosol index (AAI). It is thus of interest to investigate possible relationships between the AAI and the cloud shadow fraction and the NO 2 TVCD. We use the TROPOMI AAI product where the AAI has "positive values of the reflectance residue between an absorbing aerosol loaded atmosphere and a clear atmosphere. Negative values are associated with an atmosphere that contains more scattering particles than a clear atmosphere" (Kooreman et al.,

Results
For both the synthetic and observational data the aim is to compare the NO 2 TVCD for cloud affected cases with the "true" NO 2 TVCD and identify potential biases. This is straightforward for the synthetic data as the "true" NO 2 TVCD is known. For the observational data the "true" NO 2 TVCD is in general not known. One approach to estimate the "true" NO 2 TVCD from the observational data is discussed in section 4.2, which also include the analysis of the observational data.

Synthetic satellite data
For the synthetic data the fully cloudy pixels were not simulated due to the computational burden, see Emde et al. (2021) for details. To further filter for the presence of clouds the NO 2 TVCDs were retrieved for pixels where the radiometric cloud fraction CF r < 0.3 for all LEO and GEO cases. The retrieval tries to correct for the presence of remaining clouds by including a standard photon path length correction based on absorption by the oxygen collision pair O 2 -O 2 or the O 2 -A band as described 10 by Yu et al. (2021). An example of the retrieved NO 2 TVCD and the bias is shown in Figs. 1a and b, while the cloud shadow index is shown in Fig. 1c.
For three LEO cases the difference between the retrieved and "true" NO 2 TVCDs versus the cloud shadow index is shown in Fig. 3. The cloud shadow impact is seen to increase as the solar zenith angle increases with the number of pixels with NO 2 TVCD differences < −20% being 0.1%, 4.% and 20.3% for solar zenith angles of 20 • (Fig. 3a), 40 • (Fig. 3b) and 60 • (Fig. 3c), 15 respectively. As the solar zenith angle increases a linear relationship appears between the NO 2 TVCD difference and the CSI.
For the cloud correction of AMF the weighted radiometric cloud fraction, CF w , is used as described by Yu et al. (2021). The effect of increasing CF w is shown in Fig. 4 where the NO 2 AMF bias (1D AMF -3D AMF) is plotted as a function of CF w for all GEO and LEO geometries. The retrieval bias is binned in CF w bins with 2% steps which gives more than 5000 pixels in each bin when the CF w < 75%. For both geometries the NO 2 AMF bias is high (median of 20%) and with a large range 20 Figure 4. Distribution of the NO2 AMF bias as a function of the weighted radiometric cloud fraction, CFw, for geostationary and low earth orbiting geometries. Results are shown for the retrieval using the O2-O2 cloud correction. Results are similar for the O2-A band based cloud correction.
(0-65%) for CF w < 1%. The bias decrease to 0% when the CF w is between 1-3%. The bias range increases slightly with inreasing CF w : about 75% of the pixels have a positive bias for large CF w for GEO geometry, while there are comparatively more pixels with a negative bias for LEO geometry. The solar zenith angle is the same for the LEO and GEO cases. However, the solar azimuth angle and the sensor viewing zenith and azimuth angles are different, thus giving different sensitivity to cloud shadows for LEO and GEO geometries. This is discussed in more detail below. More than 15% of the pixels have a bias larger 5 than 20% for CF w > 50% for both GEO and LEO geometries. As discussed by Yu et al. (2021), the NO 2 retrieval error due to the use of a simple cloud correction scheme is generally less than 20% for nearly cloud-free cases (CF w < 50%). Therefore, pixels with a NO 2 retrieval bias >20% are likely to be affected by cloud shadow effects. Figure 5 shows the number of pixels that satisfy these conditions for all the LEO and GEO geometries. For the geostationary geometry, the number of pixels with retrieval bias > 20% is up to about 1000 or 10.6% (out of 9400 pixels for a single case). This number increases with surface 10 albedo and SZA, the difference between surface albedos of 0.05 and 0.2 is about 150 to 400 pixels. This number also strongly depends on the solar azimuth angle (SAA), and the difference between SAAs of 315 • and 45 • is more than 25%. For low earth orbiting geometry, the number of pixels with bias > 20% increases with surface albedo and SZA as well, and it is up to 1600 (17%) for high surface albedo of 0.2 and SZA of 60 • .
To identify the localisation of the pixels with NO 2 AMF bias > 20%, maps of the number of cases with NO 2 AMF bias 15 > 20% were made as shown in Fig. 6. Generally, the majority of pixels with large NO 2 AMF bias are found in cloud free regions close to cloud edges. Figure 7 shows maps of the maximum retrieval bias for each pixel. As above, cloudy data (CF w > 50%) are excluded from the analysis. The maximum bias is usually less than 60% over the northwest region (x=30-50, y=70-100) where there is an ice cloud with small optical thickness (< 10). In the center of the map, the bias is often higher than 100%,    the pixels around are covered by a convective cloud with large vertical extent and optical thickness > 100. Fig. 8 shows under which geometry the maximum bias is observed. Pixels with maximum biases < 20% are not shown. In general, the maximum bias is obtained at high SZA (60 • ), and for GEO cases, the bias also depends on the solar azimuth angle. Maximum biases are found on the east/west side of the cloud when the sun is in the west/east. Note that this dependency for GEO geometry is particular for the region, and thus solar and viewing conditions, studied. This dependency is not seen for the LEO geometry as 5 for LEO geometry the local daily revisit time is the same and thus also the solar azimuth.
Generally, as shown above, the NO 2 AMF retrieval bias is mostly within 20%, and the bias increases slightly with the weighted radiometric cloud fraction. The NO 2 profile and surface albedo input gives differences between the LIDORT and MYSTIC RTM on the order of 1% (Emde et al., 2021;Yu et al., 2021). Therefore, the retrieval bias is mainly due to the simple cloud correction and 3D cloud effects. For the clear pixels (CF w < 1%), there is a significant number of pixels with the 10 retrieval bias > 20% and this is due to cloud shadow effects.
To summarize the overall differences between the retrieved and "true" NO 2 TVCDs for all LEO and GEO cases, various statistics were calculated. These include the median, mean, standard deviation and skewness of the difference distribution. Furthermore, the number n of pixels with differences d: d < −20; −20 < d < −10; −10 < d < −5; −5 < d < 5; 5 < d < 10; 10 < d < 20; d > 20% were calculated. These results are summarized in Fig. 9 for the low-earth-orbiting geometries and in Fig. 10 15 for the geostationary geometries. For the low-earth-orbiting geometry the bias varies with solar and viewing zenith angles (red solid line, Fig. 9b). The average median bias is -0.5%. Between 53% (large solar zenith angle) and 89% (small solar zenith angle) of the retrieved NO 2 TVCDs are within ±10% of the true column. The number of pixels with differences < −20% Figure 9. Summary of NO2-retrieval differences for the low-earth-orbiting geometry. (a) The number of pixels with NO2 TVCD differences larger than 20% (magenta solid line) and less then -20% (blue solid line). Also shown are pixels in other difference intervals as indicated in the legend. Note that the % is with respect to the total number of pixels for which a NO2 retrieval has been done. This number is smaller than the number of pixels in the scene as cloud contaminated pixels are excluded from the retrieval. increase from about 0.1% for low albedo and high sun to up to 26% for high albedo and low sun (solid blue line, Fig. 9a).
Within each solar zenith angle interval the solar azimuth angle takes two values and the results are similar for these two angles.
The largest differences within each solar zenith angle interval are found when the viewing zenith angle is large (60 • ). The viewing azimuth angle has less influence.
For the geostationary geometry, there is an overall negative bias of about -0.9% (red line, Fig. 10b). Between 61% (large 5 solar zenith angle) and 93% (small solar zenith angle) of the retrieved NO 2 TVCD are within 10% of the "true" column. The number of pixels with differences < −20% increase from about 0.2% for low albedo and high sun to up to 22% for high albedo and low sun (Fig. 10b). The variations within each solar zenith angle interval is due to different solar azimuth angles. The skewness is always negative, indicating a skewed distribution.
These results show that the solar zenith angle is of prime importance for 3D cloud impacts and that the impact increases 10 with increasing solar zenith angle. For similar reasons, the viewing zenith angle is also of importance. Both under-and overestimates occur in pixels close to clouds. Cloud shadows are a cloud feature metric that may be used to identify affected pixels, Fig. 3. However, while for large solar zenith angles pixels affected by cloud shadows are mostly underestimated, overestimates 13 https://doi.org/10.5194/amt-2021-331 Preprint. Discussion started: 3 November 2021 c Author(s) 2021. CC BY 4.0 License. Figure 10. As in Fig. 9, but for geostationary orbit geometry. also occurr. Thus, cloud features in neighbouring pixels, such as cloud top altitude and cloud optical thickness, are also of importance.
According to box cloud simulations presented by Emde et al. (2021) and Yu et al. (2021), the cloud enhancement effect is of similar magnitude as the cloud shadow effect. For the synthetic data the enhancement effect is present, but is smaller in magnitude than the cloud shadow effect, see Figs. 9 and 10. This is most likely due to cloud enhancement affected pixels being 5 identified as cloudy and thus not analysed. Also the synthetic data indicate that the smaller the solar zenith angle is, the larger are the chances for cloud enhancements between 5-10% for low earth orbit satellites (θ = 20 − 30 • , purple dotted line Fig. 9a).
The differences are largest for large satellite viewing angles (60 • ) and may thus indicate enhancement due to the satellite partly viewing sun-illuminated cloud sides. For geostationary satellite this effect is not present (Fig. 10a). This is due to differences in the solar azimuth angles for the two geometries. However, in magnitude, this effect has a much smaller impact than cloud 10 shadows.
For comparison Lorente et al. (2017) estimated the structural uncertainty (differences in retrieval methodology) in the tropospheric NO 2 AMF to be 42% over polluted regions and 31% over unpolluted regions, and mostly driven by the uncertainty in the a-priori NO 2 profile, cloud properties and surface albedo. Furthermore, Yu et al. (2021) showed that the NO 2 AMF bias using the simple cloud correction is generally within 20% for cases with CF w < 50%. Yu et al. (2021) found that the retrieval 15 bias significantly increases for high albedo (albedo=0.3) for 1D cloud layer simulations, whereas for box-cloud sensitivity studies, the 3D bias decreases with increasing surface albedo. The bias of the NO 2 retrieval has two origins: the simplistic treatment of clouds in the retrieval algorithm, and 3D effects. Due to the different results for the 1D and box-cloud results, the high bias for high albedo cases is believed to be due to the simplistic treatment of clouds in the retrieval algorithm. The

Cloud shadow band cases
For the observational data the "true" NO 2 TVCD unaffected by clouds is not known. In order to have a reference or "true" NO 2 TVCD to compare with, we look at neighbour pixels in a 3×3 pixel matrix where the pixel of interest is in the centre. The "true" NO 2 TVCD is then taken to be the average of cloudfree neighbours with NO 2 retrieval quality value > 0.95. Obviously this choice of "true" NO 2 TVCD has its problems including that neighbours may also be affected by clouds and that NO 2 10 TVCD may have large spatially gradients, see for example Fig. 2b which shows the NO 2 TVCD, and Fig. 2c which shows the percentage difference of the NO 2 TVCD to the area median NO 2 TVCD.
Based on the findings from the analysis of the synthetic data we searched TROPOMI and VIIRS data for cases with cloud shadows and large solar zenith angles. An example of a cloud shadow band is shown in Fig. 11. The NO 2 TVCD, Fig. 11b, is markedly higher in the Amsterdam area, but is otherwise overall slowly varying over the region. The H-metric, Fig. 11c, is 15 sensitive to inhomogeneous clouds, but also to variations in the surface albedo, compare Fig. 11a and 11c. The geometric cloud fraction, based on the VIIRS cloud mask, shows larger variability than the average radiometric cloud fraction CF V IIRS r , as expected, compare Fig. 11d and 11e and see also the distributions of these quantities in Fig. 11h. The cloud shadow fraction, Fig. 11f, is between 0.2 and 0.5 in the area with scattered clouds. The cloud shadow fraction also clearly delineates the large cloud structures. The NO 2 TVCD versus the cloud shadow fraction is shown in Fig. 11g. For a cloud shadow fraction < 0.5 20 the NO 2 TVCD decrease with increasing cloud shadow fraction. However, the decrease is well within the variability of the NO 2 TVCD (the red line with error bars indicates the average NO 2 TVCD ±standard deviation). The distribution of the Hmetric for various CF V IIRS r is shown in Fig. 11i. For CF V IIRS r between 0.4 and 0.6 (green bars) maxima are found both for low (scattered clouds) and high (homogeneous clouds or surface) H-metrics. Thus indicating that the CF V IIRS r may not unambigiously be used to identify scattered cloud cases. For the absorbing aerosol index no obvious dependence between the 25 AAI and the NO 2 TVCD is present for this region, see Appendix A for further details.
For the cloud shadow band pointed to by the red arrow in Fig. 11a, further analysis is provided in Fig. 12. A RGB zoom in of the cloud shadow band is shown in Fig. 12b where red marks indicate pixels with cloud shadow. The VIIRS cloud top height and the TROPOMI NO 2 TVCD are shown in Figs. 12c and 12d. The cloud top height, cloud optical thickness and the cloud fraction is shown for row 267 in Fig. 12e. The cloud shadow is just north of a cloud band with optical thickness up 30 to 10 (slant optical thickness of about 15) and an altitude between 9-10 km. North of the cloud shadow band there are also some clouds at the same height but with a smaller optical thickness of around 1.5-2. These thinner clouds are not easily seen in the VIIRS RGB, Fig. 11a, but are present in the VIIRS cloud mask and thus the geometric cloud and averaged radiometric 262-269 and the average of these. A shadow pixel is defined as having a CSF > 50% and CF w < 50%. For row 263 no pixels satisfied this criteria and does no data is shown in the shadow region for this row. Pixels to be included in the regions north and south of the shadow band (up to 4 TROPOMI pixels north and south of shadow band), where required to have CSF < 25% and CF w < 50%. For the cloud shadow band the NO 2 TVCD is on average reduced by 17%. There is no clear dependence of the NO 2 difference on the cloud shadow fraction. The cloud optical thickness of the cloud causing the cloud shadow varies little, 5 the average cloud optical thickness being 9.0±1.6. Thus it is not possible to say anything about the dependence of the NO 2 difference on the cloud optical thickness for this case.
Another example of a cloud shadow band and 3D cloud effects on NO 2 retrieval is shown in Fig. 13. The VIIRS RGB image ( Fig. 13a) shows the cloud coverage over Northern Germany on December 30, 2019, when the solar zenith angle during the VIIRS and TROPOMI overpass time is larger than 70 • . The cloud shadow is just north of a large cloud band, and most of the 10 northern region is completely cloudless, which is similar to the idealized box cloud cases presented by Emde et al. (2021) and Yu et al. (2021). The cloud height from west to east is 1 km to 8 km (Fig. 13c). The cloud optical thickness is not available for this case, but from the RGB the cloud is clearly optically thick enough to make the ground not visible from space. To look for cloud shadow effects, we select the region within the red box in Fig. 13a). In Fig. 13e is shown the VIIRS reflectance and cloud top height for TROPOMI row 394 as a function of latitude. The pixels are identified as cloudy, cloud shadow and clear 15 based on the reflectance; with the reflectance being about 0.25 over the clear region, down to 0.18-0.24 in the cloud shadow, and higher than 0.25 for cloudy pixels. There are four TROPOMI pixels in the cloud shadow where the NO 2 TVCD is low by 20-60% compared with the NO 2 TVCD to the north and south of the cloud shadow (cloud and clear pixels). There is a slight difference in the NO 2 retrieval using different cloud corrections, Fig. 13f. In order to reduce the influence of the signal to noise ratio, the NO 2 TVCD is averaged over cloudless (four pixels near the cloud shadow), cloud shadow and cloudy pixels (four 20 pixel near the cloud shadow) for rows 393 to 398, as shown in Fig. 13g. For the cloudy pixels, the difference is within 10% between NO 2 averaged over all the pixels and NO 2 averaged for the pixels with high quality retrieval (CF w < 50%). All cases show that the NO 2 TVCD in the cloud shadow is lower by 8-46% (average of 25%) compared with the NO 2 TVCD around the shadow. Here, the NO 2 TVCD around the shadow represents NO 2 retrieval unaffected by 3D cloud, which is an average of NO 2 over cloudy and cloudless pixels, in order to reduce the impact of spatial variation of NO 2 .

25
If it is assumed that the clouds are the main reason for the variations in the NO 2 TVCD over the cloud shadow bands, then these cases are examples of how cloud shadows give underestimates of NO 2 TVCD, in agreement with the theoretical idealized box cloud results presented by Emde et al. (2021) and Yu et al. (2021).

Cloud effects for general cases
As shown in Figs. 9 and 10 for the synthetic data, the NO 2 bias may be on the order of tens of percent and are largest for large 30 solar zenith angles. Hence, for the months of October 2018 and March 2019 when the solar zenith angle is between 50-60 • for the study region, we looked for NO 2 TVCD biases due to clouds. For October 2018 and March 2019 this amounted to a total of 1,023,081 TROPOMI pixels. To quantify possible cloud shadow effects and factors that impact it, TROPOMI pixels with NO 2 retrieval data quality value >0.95 and cloud shadows were selected for further analysis. A NO 2 retrieval with the data quality value >0.95 was reported for 367,584 (36%) of the pixels. Of these, 129,180 (35%) were affected by cloud shadows according to the VIIRS cloud shadow product. Note that this number pertains to months of the year where we expect cloud shadow effects to be large. For months of the year with smaller solar zenith angles this number may well be smaller. The pixels affected by cloud shadows were further analysed to understand which parameters that have the largest impact on the NO 2 TVCD. In order to have a reference NO 2 TVCD to compare the centre pixel value with, only centre pixels which have one or more cloudfree neighbours with NO 2 retrieval quality value > 0.95, were included. This restriction reduced the number of TROPOMI pixels to be analysed from 129,180 to 39,011. For these pixels the difference between the centre pixel NO 2 TVCD and the average of the NO 2 TVCD in the cloud free neighbours (∆NO 2 ) is shown as a function of cloud height maximum of neighbour pixels in Fig. 14. For this data subset clouds are mostly found at high (about 10,000 m) and low (about 1,900 m) altitudes, with relatively more clouds at the higher altitudes. This is qualitatively in agreement with Noel et al. (2018). They reported cloud fractions as 5 estimated from a lidar on board the International Space Station. For Europe (their Fig. 5b) they report a similar cloud vertical behaviour, albeit for a different time of year (July, June and August). The subset of TROPOMI pixels presented in Fig. 14 shows both high and low NO 2 TVCD biases with a median NO 2 TVCD bias of -1.4% and maximum of the Gaussian kernel density probability density function estimate maximum of -0.7%. Thus, for this subset of TROPOMI pixels no signifcant cloud shadow effect was visible in the NO 2 TVCD. As stated above, no true NO 2 TVCD is available as for the synthetic data and 10 this may be the reason for the lack of a clear cloud shadow effect in this data subset.
This subset of data were further binned according to the maximum cloud height of neighbour pixels and the maximum slant cloud optical thickness (SCOT) of neighbour pixels as given in Table 1. The maximum slant cloud optical thickness is the Maximum slant cloud optical thickness of neighbour pixels (SCOT) 0, 3,5,7,10,15,300 neighbour pixel maximum cloud optical thickness divided by the cosine of the solar zenith angle for that pixel. Furthermore, only TROPOMI pixels where the maximum slant cloud optical thickness and maximum cloud top height are in the same neighbour pixel were included. This reduced the data set to 18,029 pixels.
The NO 2 TVCD bias density as a function of the cloud shadow fraction for the SCOT bins and maximum cloud heights is shown in Fig. 15. As the cloud height increases (from bottom row to top row in Fig. 15) the cloud shadow fraction increase 5 because generally the cloud shadow within a pixel geometrically increase with cloud height when the cloud is larger than the pixel size. For the low clouds the median NO 2 TVCD bias varies between -1.2 to 0.2% for SCOT<15. For SCOT>15 the median NO 2 TVCD bias is -2.9%. For the medium height clouds the median NO 2 TVCD bias for SCOT>15 increase in magnitude to -5.8%. For smaller SCOT the median NO 2 TVCD bias is negative and varies between -0.8 and -3.4%. For the high clouds the median NO 2 TVCD bias is between -0.3 and 1.6% for SCOT<7. For larger SCOT the median NO 2 TVCD 10 bias is -0.9, -8.7 and -15.1% (7<SCOT<10, 10<SCOT<15, 15<SCOT<300).
For all cloud heights the median NO 2 TVCD bias is negligible for SCOT<7. For low and medium clouds there is a negative bias for the largest SCOT. For the high clouds the bias is pronounced for SCOT>10. Thus, the median NO 2 TVCD bias increases with cloud height and with slant cloud optical thickness. It is noted that for individual pixels the bias may be larger and both positive and negative. We use a neighbour pixel as the reference NO 2 -value. This assumes that only the cloud shadow 15 effect is the reason for the NO 2 TVCD bias. In reality there are horizontal gradients in the NO 2 TVCD for numerous other causes as well, including local emissions and wind transport differences.
Finally it is noted that for the 129,180 pixels which were affected by cloud shadows, the number of pixels with cloud free neighbours where | ∆NO 2 |> 20% is 45.7%. Under the assumption that this number is applicable to all TROPOMI pixels affected by cloud shadows, about 16% of TROPOMI pixels for which NO 2 TVCD retrievals with a high quality value were 20 made, may be impacted by cloud effects larger than 20% for solar zenith angles between 50-60 • . Of these about half of the NO 2 TVCD are overestimated and half underestimated. The underestimate is clearly linked to cloud shadow effects. The overestimate may be due to in-scattering or horizontal gradients in NO 2 concentrations and thus wrong cloud shadow free NO 2 TVCD reference. These two processes may, however, not be distinguished in the observed data set.

25
In this study we have investigated the impact of 3D clouds on NO 2 TVCD retrievals from UV-VIS sounders. Both synthetic and observational data have been used to identify and quantify possible biases in NO 2 TVCD retrievals. The synthetic data were based on high-resolution LES results which are input to the MYSTIC 3D radiative transfer model. The simulated visible spectra  for low-earth orbiting and geostationary geometries were analysed with standard NO 2 retrieval methods including operational cloud corrections. Profiles of NO 2 for polluted conditions were considered as cloud shadow effects are not important for background NO 2 conditions. For the observational data the NO 2 products from TROPOMI were used while VIIRS provided high spatial resolution cloud data. Both single cases and overall statistics were calculated. The main findings are: -The following metrics were identified as being the most important in identifying 3D cloud impacts on NO 2 retrievals: 5 cloud shadow fraction, cloud top height, cloud optical thickness, solar zenith and viewing angles.
-Analysis of the synthetic data show that for low-earth and geostationary orbit geometries, 89 and 93%, respectively, of the retrieved NO 2 TVCD are within 10% of the actual column for small solar zenith angles. For large solar zenith angles the numbers decrease to 53 and 61%.
-The synthetic data shows that in general the NO 2 TVCD bias is slightly larger for low-earth orbiting geometries than 10 for geostationary geometries. This is due to differences in viewing geometry where, for the mid-latitude targets studied here, the sun-target-detector geometry overall sees fewer cloud shadows for geostationary geometry.
-For a solar zenith angle less than about 40 • the synthetic data show that the NO 2 TVCD bias is below 10%. For larger solar zenith angles both synthetic and observational data show NO 2 TVCD bias on the order of tens of %.
-For clearly identified cloud shadow bands in the observational data, the NO 2 TVCD appears low-biased when the cloud 15 shadow fraction > 0.0 compared to when the cloud shadow fraction is zero. If it is assumed that the clouds are the main reason for the variations in the NO 2 TVCD over the cloud shadow band, i.e. the NO 2 field is assumed to be horizontally homogeneous, then these cloud shadow band cases are examples of how cloud shadows give underestimates of NO 2 TVCD, in agreement with the theoretical findings.
-For solar zenith angles between 50-60 • , about 16% of TROPOMI pixels (spatial resolution 3.5×7 km 2 at nadir) for 20 which NO 2 retrievals with a high quality value were made, were found to be impacted by cloud effects larger than 20%.
The above conclusions suggest that further work is required on 3D cloud radiative impacts. For this future work the following topics may be considered: -As shown in this study, 3D radiative simulation of high-resolution satellite instrument spectra with realistic 3D cloud input is fully achievable. However, to cover all possible atmospheric situations on Earth, possibly more high-spatial 25 resolution LES results may be needed. Furthermore, 3D RT modelling is demanding on computer resources and requires careful interpretation and analysis of the output. It is expected that such modelling efforts may prove useful not only for trace gas retrieval algorithm studies, but also for cloud detection and cloud microphysics retrievals and more.
- shadow product (large changes between versions) suggests that independent verification with ground measurements may be of use. Such validation is non-trivial and possibly require new experimental approaches for measurements of both cloud shape and trace gas spatial variation at sub-pixel resolution.
Code availability. The libRadtran software used for the radiative transfer simulations is available from www.libradtran.org. The QDOAS software for DOAS retrieval of trace gases is available from uv-vis.aeronomie.be/software/QDOAS.

Appendix A: Absorbing aerosol index
Recently Kooreman et al. (2020) identified and explained AAI large scale effects such as the cloud bow, sunglint and viewing zenith angle dependence. In addition they reported small-scale negative AAI values in partly cloudy areas and attributed this 10 to 3D cloud structures casting shadow, but left the in-depth analysis for a future study, As cloud shadow impact NO 2 TVCD retrievals it is of interest to investigate possible relationships between the AAI and the cloud shadow fraction and the NO 2 TVCD. For the partly cloudy scene in Fig. 11, the TROPOMI AAI from 380 and 340 nm is shown in Fig. A1a. Overall the AAI is negative indicating the absence of absorbing aerosol, but the presence of scattering particles. Indeed, over land the AAI is more negative over cloudy pixels, compare Fig. 11d