The importance of surface reflectance anisotropy for cloud and NO₂ retrievals from GOME-2 and OMI

The angular distribution of the light reflected by the Earth's surface influences top-of-atmosphere (TOA) reflectance values. This surface reflectance anisotropy has implications for UV/Vis satellite retrievals of albedo, clouds, and trace gases such as nitrogen dioxide (NO₂). These retrievals routinely assume the surface to reflect light isotropically. Here we show that cloud fractions retrieved from GOME-2A and OMI with the FRESCO and OMCLDO2 algorithms have an east–west bias of 10 % to 50 %, which are highest over vegetation and forested areas, and that this bias originates from the assumption of isotropic surface reflection. To interpret the across-track bias with the DAK radiative transfer model, we implement the bidirectional reflectance distribution function (BRDF) from the Ross–Li semi-empirical model. Testing our implementation against state-of-the-art RTMs LIDORT and SCIATRAN, we find that simulated TOA reflectance generally agrees to within 1 %. We replace the assumption of isotropic surface reflection in the equations used to retrieve cloud fractions over forested scenes with scattering kernels and corresponding BRDF parameters from a daily high-resolution database derived from 16 years' worth of MODIS measurements. By doing this, the east–west bias in the simulated cloud fractions largely vanishes. We conclude that across-track biases in cloud fractions can be explained by cloud algorithms that do not adequately account for the effects of surface reflectance anisotropy. The implications for NO₂ air mass factor (AMF) calculations are substantial. Under moderately polluted NO₂ and backward-scattering conditions, clear-sky AMFs are up to 20 % higher and cloud radiance fractions up to 40 % lower if surface anisotropic reflection is accounted for. The combined effect of these changes is that NO₂ total AMFs increase by up to 30 % for backward-scattering geometries (and decrease by up to 35 % for forward-scattering geometries), which is stronger than the effect of either contribution alone. In an unpolluted troposphere, surface BRDF effects on cloud fraction counteract (and largely cancel) the effect on the clear-sky AMF. Our results emphasise that surface reflectance anisotropy needs to be taken into account in a coherent manner for more realistic and accurate retrievals of clouds and NO₂ from UV/Vis satellite sensors. These improvements will be beneficial for current sensors, in particular for the recently launched TROPOMI instrument with a high spatial resolution.

1 Introduction

Nitrogen dioxide (NO₂) in the lower troposphere is an important constituent of air pollution. In Europe, the annual mean NO₂ concentration limit value (40 µg m⁻³) is still widely exceeded, exposing 30 million people to poor air quality with known harmful health effects (EEA, 2016). In combination with other pollutants and sunlight, chemical and physical transformations of nitrogen oxides ($\mathrm{NO}+{\mathrm{NO}}_{\mathrm{2}}={\mathrm{NO}}_x$) lead to the formation of particulate matter and ozone smog, further impacting public health, ecosystems, and climate. Satellite measurements of tropospheric NO₂ column densities provide much better spatial coverage than ground-based sensors, and they have been used to monitor trends and to estimate NO_x emissions and NO₂ surface concentrations (e.g. Richter et al., 2005; Martin et al., 2003; Lamsal et al., 2008). The spatial resolution of the satellite instruments and their retrievals is improving such that the observed NO₂ pollution can now be traced back to emissions from individual cities, power plants, and transportation sectors. The uncertainty of satellite NO₂ retrievals is considerable and mainly related to the adequacy of the assumptions made on the state of the atmosphere. We recently estimated the structural uncertainty from an ensemble of NO₂ retrievals to be of the order of 30 %–40 % (Lorente et al., 2017). An important component of this uncertainty is how surface properties (usually from an external database) are taken into account, and how errors in the external database propagate in the air mass factor (AMF) calculations. This is not straightforward, because the AMF calculation directly depends on surface properties under clear-sky circumstances, and indirectly via cloud parameters retrieved for the same scene by the cloud algorithm.

Surfaces reflect light differently in each direction, and the angular distribution of the reflected light influences top-of-atmosphere (TOA) reflectance levels measured by satellite instruments that monitor atmospheric composition. Therefore, surface reflectance anisotropy influences retrievals of surface albedo, trace gases, aerosols, and clouds from satellite instruments like the Global Ozone Monitoring Experiment 2 (GOME-2) and the Ozone Monitoring Instrument (OMI). In surface albedo, cloud, aerosols, and trace gas retrievals, the surface is often assumed to be Lambertian: an idealised surface that reflects light isotropically (e.g. Kleipool et al., 2008; Veefkind et al., 2016; Torres et al., 2007; Boersma et al., 2011). This assumption implies that the geometry-dependent scattering properties of the reflecting surface are ignored.

The Lambertian-equivalent reflectivity (LER) climatologies represent the albedo of the Lambertian surface in the radiative transfer simulations for cloud retrievals and trace gas retrievals. In constructing such climatologies (e.g. monthly climatologies), a large ensemble of measurements taken over a scene over multiple years is analysed statistically, and based on the lower 1 % percentile reflectance, an inversion is done to retrieve the surface reflectance (e.g. Koelemeijer et al., 2003; Kleipool et al., 2008; Tilstra et al., 2017). Depending on the exact viewing and illumination geometry, however, the surface may appear darker or brighter. Taking the lower 1 % percentile reflectances therefore skews the distribution of retrieved albedo values to those scenes that appear darker from space. Using these climatologies therefore fails to represent any surface reflectance anisotropic effects on TOA reflectance simulations for the widely varying subset of viewing geometries encountered along a satellite orbit. We may expect cloud retrievals to be directly affected by the assumption of a Lambertian surface in the radiative transfer: if a scene is brighter than predicted by the biased climatology, cloud fractions will be overestimated. Trace gas retrievals are affected by the Lambertian assumption in the calculation of the AMF: directly via the clear-sky AMF and also indirectly because the retrieved cloud parameters are used to correct for the possible presence of residual clouds via the independent pixel approximation (IPA) method. In the IPA, the cloud radiance fraction weighs the clear-sky and cloudy parts of the scene for the calculation of the overall AMF and vertical column density (VCD) (e.g. Martin et al., 2002).

The angular distribution of the reflected light by a surface is represented mathematically by the bidirectional reflectance distribution function (BRDF) (Nicodemus et al., 1992). Anisotropy is a fundamental physical property of surface reflectance, so in order to fully represent the geometry-dependent surface scattering properties in cloud and trace gas retrievals, surface BRDF has to replace the isotropic Lambertian albedo. Some studies have already shown that surface BRDF effects are important for NO₂ and cloud retrievals. Zhou et al. (2010) found that, after including surface BRDF over Europe, differences in NO₂ columns were higher in November (20 %) compared to July (3 %), when the solar zenith angle is small. Noguchi et al. (2014) studied surface BRDF effects on the diurnal cycle of clear-sky geostationary measurements over Japan, and found that whether NO₂ columns are under- (−15 %) or overestimated (+9 %) depends on the specific geometry of the measurement. To also address the need to include surface BRDF effects on cloud algorithms, Lin et al. (2014, 2015) updated the POMINO retrieval over China. Changes in surface reflectance led to changes in opposite sign in cloud fraction (±0.05) and more complex effects on cloud pressure, with an overall change of 10 % in NO₂ columns that were regionally and seasonally dependent. Vasilkov et al. (2017) created a geometry-dependent surface LER (GLER) product and applied it to OMI NO₂ and cloud retrievals. They found relatively small effects on retrieved cloud parameters and relatively high differences in retrieved NO₂ columns (up to 50 %) driven by GLER values, being on average smaller than the original LER.

These studies show that surface BRDF effects depend on the specific geometry (hence local time and season) and surface characteristics of the individual measurements and that averaging over many pixels results in smaller differences. They analysed mostly clear-sky scenes (i.e. no clouds present) or scenes with very low cloud fractions (i.e. lower than 0.2), and they did not consider how surface reflectance anisotropy affects the radiative transfer in the atmosphere and TOA reflectances. The indirect effect of biased cloud parameters on NO₂ retrievals combined with effects on clear-sky AMFs received less attention. Vasilkov et al. (2017) and Lin et al. (2015) addressed the indirect effects and showed that the effects on cloud parameters could enhance or compensate for the direct effect on clear-sky AMFs.

Here we study the effect of surface reflectance anisotropy on surface LER climatologies, on cloud fraction retrievals and on NO₂ tropospheric AMFs, covering the essential steps from the TOA reflectance to retrieve the final NO₂ tropospheric column in a partly cloud-covered pixel. We present observational evidence that surface reflectance anisotropy skews LER climatologies to the darkest scenes corresponding with forward-scattering geometries (Sect. 2). We demonstrate that using these LER climatologies in cloud retrieval algorithms leads to considerable across-track biases in cloud fractions, especially for the O₂-A band (GOME-2) but also for the 477 nm O₂–O₂ band (OMI). In Sect. 3 we describe our extension of the DAK radiative transfer model (RTM) to include surface reflectance anisotropy with the Ross–Li BRDF semi-empirical model. We validate DAK with state-of-the-art radiative transfer models and evaluate TOA reflectance simulations including surface BRDF at the relevant wavelengths for cloud and NO₂ retrievals. In Sect. 4 we study the consequences of including surface BRDF in the calculation of effective cloud fractions. In Sect. 5 we study the surface BRDF effects on NO₂ AMFs, and how these, in combination with the effects on cloud fractions, affect tropospheric NO₂ retrievals in cloudy scenes. We end with conclusions and an outlook.

2 Evidence of the influence of surface reflectance anisotropy on LER climatologies and cloud retrievals

Surface LER climatologies are commonly used as boundary conditions for cloud and trace gas retrievals, (e.g. De Smedt et al., 2018; Bucsela et al., 2013). Earlier instruments like GOME and SCIAMACHY have very coarse pixels (40×320 km² and 60×30 km² respectively) and a narrow swath (960 km). For retrievals from these instruments, the Lambertian assumption can be justified as surface BRDF effects are likely to smooth out over the large and heterogeneous pixels. Newer instruments like OMI, and especially TROPOMI, have higher spatial resolution (13×24 km² and 3.5×7 km² respectively) and a wider swath (up to 2600 km), i.e. a wider range of viewing directions, and therefore the surface BRDF effects become more relevant. One of the advantages of using LER climatologies is that they have been derived from measurements of the satellite instrument itself, for example, Koelemeijer et al. (2003) climatology for GOME and SCIAMACHY, Kleipool et al. (2008) for OMI, and Tilstra et al. (2017) for GOME-2 and SCIAMACHY. In constructing these climatologies, it is assumed that the surface reflectance is fully isotropic, and the angular dependence of reflected light is neglected. In the minimum-LER method, surface LER values are the 1 % cumulative values retrieved from the histogram of Earth reflectance over a specific scene in a climatological period. The presence of clouds increases TOA reflectance compared to clear-sky scenes; therefore taking the minimum LER ensures that the surface LER values represent mostly cloud-free scenes. However, if for particular viewing and illumination geometries, the TOA reflectance appears brighter, then those measurements will not be included in the climatology. This will introduce an intrinsic (but not explicit) viewing-angle dependency in the surface LER climatologies derived from the satellites. That measurements taken under different viewing geometry configurations are used to create the climatologies does not mean that climatological surface LER values are representative for typical or average-viewing geometries of the instruments.

Figure 1 (a, b) Map of GOME-2A minimum surface LER climatology (2007–2013) for March (a) at 494 nm and (b) at 772 nm over Amazonia (lat: 5^∘ N–10^∘ S, long: 60–70^∘ W, upper panel) at $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ resolution. (c, d) Directional dependence of surface LER climatology (2007–2013) derived from individual measurements along the swath: east (E) for the eight easternmost pixels, nadir (N) for the eight centre pixels and west (W) for the eight westernmost pixels. The horizontal line represents the surface LER using the full swath which is the average of the surface LER in panels (a) and (b).

Figure 1a, b show the minimum surface LER climatology (2007–2013) at 494 and 772 nm for March over Amazonia for GOME-2 on-board MetOp-A (hereafter GOME-2A) from Tilstra et al. (2017). This climatology was derived from the 1 % cumulative reflectances gathered irrespective of viewing geometry. Figure 1c, d show the directional dependence of the minimum surface LER climatologies derived using the same GOME-2A measurements but now discriminating between different geometries. For this purpose, the 24 measurements in 24 different viewing directions along the GOME-2 swath are considered independently, and then the minimum LER method is applied for each viewing direction as in the derivation of the original full swath climatology. We consider east as the eight most eastward measurements, nadir as the centred ones and west as the eight most westward measurements. The average relative azimuth angle (RAA) characterises the light-scattering regime for measurements in each region: low RAA in the east (RAA = 16^∘) corresponds to forward scattering and high RAA (RAA = 164^∘) in the west corresponds to backward scattering.

In Fig. 1c, d the horizontal line is the original surface LER value obtained using all GOME-2A measurements over Amazonia, and corresponds to the average value over the box in Fig. 1a, b ($\stackrel{\mathrm{‾}}{A_{\text{LER}}}=\mathrm{0.028}$ at 494 nm and $\stackrel{\mathrm{‾}}{A_{\text{LER}}}=\mathrm{0.21}$ at 758 nm). The albedo is lower at 494 nm because of the absorption of light by chlorophyll. At 758 nm, if only measurements of the eastern part (E) of the orbit are used to construct the climatology, the average surface LER value is slightly lower ($\stackrel{\mathrm{‾}}{A_{\text{LER}}}=\mathrm{0.19}$) than the value using all measurements. If only nadir (N) measurements are used, the average surface LER increases ($\stackrel{\mathrm{‾}}{A_{\text{LER}}}=\mathrm{0.24}$), and for the westernmost measurements (W), the average surface LER increases to almost twice the original value ($\stackrel{\mathrm{‾}}{A_{\text{LER}}}=\mathrm{0.37}$). Using the full swath climatology thus implies a slight overestimation of the surface reflectance for the eastern measurements and a strong underestimation for nadir and western measurements. This systematic effect is a consequence of the directional signature of the surface reflectance. In the backward-scattering direction, canopy surfaces generally appear brighter than in the forward-scattering direction (e.g. Camacho-de Coca et al., 2004). This effect is strongest in the near-infrared (NIR, 0.7–2.5 µm) spectral range, where the atmosphere is more transparent, and over vegetation, which has non-isotropic elements (e.g. dense trees with heterogeneous leave orientation and shadowing effects). The effect in the surface LER climatologies also appears over non-vegetated regions and at shorter wavelengths. However, it is not as strong because stronger Rayleigh scattering tends to smooth out the sensitivity to surface effects and because land surfaces are darker at shorter wavelengths. Non-vegetated areas are usually more isotropic than vegetated areas. Because these biased surface LER climatologies are used in cloud retrievals (e.g. FRESCO Wang et al., 2008, O₂–O₂ Veefkind et al., 2016), we anticipate a substantial effect on the retrieved cloud properties and as a consequence on trace gas column retrievals that use cloud parameters retrieved in the NIR, where sensitivity to surface anisotropy is strong (such as GOME-2).

Figure 2 Box plot of cloud fractions retrieved with (a) FRESCO cloud retrieval for GOME-2A for March 2008 and (b) OMCLDO2 cloud retrieval for OMI for September 2005 for east–nadir–west measurements over Amazonia (lat: 5^∘ N–10^∘ S, long: 60–70^∘ W). Black triangles correspond to mean values, red lines to median. The box represents 25th and 75th percentiles and the dashed lines the minimum and maximum values.

Indeed, we find that cloud fractions retrieved with FRESCO cloud retrieval from GOME-2 measurements are affected by the across-track bias in the surface LER climatology. FRESCO retrieves cloud properties in the O₂ A-band near 760 nm. Figure 2a shows that over Amazonia (in March 2008) FRESCO cloud fractions are generally lower for the eastern measurements than for nadir and western measurements. This dependency can be explained by the directional biases in the surface LER (Fig. 1d). In the nadir and west measurements, the surface LER is underestimated and the retrieval compensates these overestimating cloud fractions in order to match observed TOA reflectances.

This results in higher mean effective cloud fractions for the nadir ($\stackrel{\mathrm{‾}}{c_{\mathrm{eff}}}=\mathrm{0.50}$) and west ($\stackrel{\mathrm{‾}}{c_{\mathrm{eff}}}=\mathrm{0.66}$) measurements compared to the east measurements ($\stackrel{\mathrm{‾}}{c_{\mathrm{eff}}}=\mathrm{0.33}$). The east–west bias ($\mathrm{100}\cdot (\stackrel{\mathrm{‾}}{c_{\mathrm{eff},\mathrm{W}}}-\stackrel{\mathrm{‾}}{c_{\mathrm{eff},\mathrm{E}}})/\stackrel{\mathrm{‾}}{c_{\mathrm{eff},\mathrm{W}}}$) in the cloud fraction depends on the time of the year and the location. It is not only present over forested areas (i.e. Amazonia (50 %), Equatorial Africa (42 %) in March 2008) but also occurs over other regions (e.g. over Europe (25 %) and Asia (10 %), not shown). Furthermore, in the ensemble of western measurements in most of the regions there are very few cloud fraction values lower than 0.2. This directly impacts trace gas retrievals, because a cloud fraction of 0.2 is often used as a threshold above which it is considered difficult to retrieve tropospheric NO₂ columns (cloud screening effect). This bias in FRESCO cloud fractions is significantly higher than the cloud fraction retrieval uncertainty estimates of 0.05 due to surface albedo uncertainty (Koelemeijer et al., 2002), which underlines the need to correct for surface BRDF effects in cloud retrievals in the O₂-A band.

The OMCLDO2 cloud retrieval from the OMI retrieves cloud properties in the O₂–O₂ band around 470 nm and uses the Kleipool et al. (2008) surface LER climatology, which is based on the same principles as the climatology used in FRESCO. Cloud fractions from the OMI instrument retrieved with the OMCLDO2 algorithm also show a west–east¹ bias (Fig. 2b) over Amazonia (September 2005) of 26 % ($\stackrel{\mathrm{‾}}{c_{\mathrm{eff},\mathrm{W}}}=\mathrm{0.31}$, $\stackrel{\mathrm{‾}}{c_{\mathrm{eff},\mathrm{E}}}=\mathrm{0.42}$) and around 15 % over other regions. The effect is weaker than for GOME-2A (50 % vs. 26 % over Amazonia) but still substantial. The bias shown here is slightly higher than the cloud fraction retrieval uncertainty estimate, which is always below 0.1 (Acarreta et al., 2004), suggesting that the bias in the O₂–O₂ retrieved cloud fractions can be significant depending on location and time of the year. Because OMI angles are larger than for GOME-2, a larger effect could be expected if the O₂-A band was also used to retrieve clouds from OMI measurements.

We have shown that surface BRDF effects result in a distinct across-track bias in surface LER climatologies and cloud fractions retrieved from satellite instruments and that the effect is highly relevant in the NIR and in the visible. Errors in cloud fraction and surface albedo are the most important sources of tropospheric AMF errors (Boersma et al., 2004), so we expect a strong impact on tropospheric NO₂ retrievals. In the following section we describe how to account for surface anisotropy in the radiative transfer model DAK. Then, we study how cloud fraction and NO₂ AMFs in the framework of cloud and trace gas retrievals are affected by the assumption of a Lambertian surface compared to a realistic anisotropically reflecting surface.

3 Reflectance simulations with surface BRDF in DAK

3.1 Definition of BRDF

The amount of radiation reflected by a surface in a certain direction depends on the direction of the incident irradiance and on the direction in which the reflected radiance is observed. The surface bidirectional reflectance distribution function (BRDF) is a function that characterises the directional reflecting properties of a surface. The surface BRDF mathematically describes the angular distribution of the surface reflectance: R_g as a function of the illumination direction (incident, θ^′, φ^′) and viewing direction (reflected, θ, φ) (see Fig. 3). It is expressed as the ratio of the reflected radiance in a certain direction (dL) and the incident irradiance from a particular direction (dE^′) (Nicodemus et al., 1992):

$$\begin{array}{}\text{(1)}& {\mathbf{R}}_{\mathrm{g}}({\mathit{\theta }}^{\prime },{\mathit{\phi }}^{\prime };\mathit{\theta },\mathit{\phi })=\mathrm{d}L({\mathit{\theta }}^{\prime },{\mathit{\phi }}^{\prime };\mathit{\theta },\mathit{\phi })/\mathrm{d}{E}^{\prime }({\mathit{\theta }}^{\prime },{\mathit{\phi }}^{\prime }).\end{array}$$

The zenith angle (θ) and the azimuth angle (φ) define the direction of incidence (${\mathit{\theta }}^{\prime },{\mathit{\phi }}^{\prime }$) and reflectance (θ,φ), as sketched in Fig. 3.

The albedo of a surface is generally defined as the ratio of the irradiance reflected by a surface area into the whole hemisphere and the irradiance incident on the surface with hemispherical angular extent (i.e. coming from all directions) (Schaepman-Strub et al., 2006). For particular illumination and viewing conditions, black-sky and white-sky albedo are defined, and they are obtained through hemispherical integration of the surface BRDF. The black-sky albedo is the albedo without diffuse component in the incident irradiance; i.e. the illumination of the surface is from a single direction. It is defined as the integral of the BRDF over the reflection hemisphere of 2π steradian (directional–hemispherical reflectance):

$$\begin{array}{}\text{(2)}& A_{\mathrm{bs}}({\mathit{\theta }}^{\prime },\mathit{\lambda })=\underset{\mathrm{0}}{\overset{{}^{\mathrm{2}\mathit{\pi }}}{\int }}\underset{\mathrm{0}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}{\mathbf{R}}_{\mathrm{g}}({\mathit{\theta }}^{\prime },{\mathit{\phi }}^{\prime };\mathit{\theta },\mathit{\phi })\cos \mathit{\theta }\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\phi }.\end{array}$$

In the particular case when there is only a diffuse isotropic component in the incident irradiance, the white-sky albedo can be defined as the integral of the surface BRDF over both the incident and reflection hemispheres (bihemispherical reflectance):

$$\begin{array}{ll}\text{(3)}& {\displaystyle }& {\displaystyle }{A}_{\mathrm{ws}}\left(\mathit{\lambda }\right)={\displaystyle }& {\displaystyle }\underset{\mathrm{0}}{\overset{{}^{\mathrm{2}\mathit{\pi }}}{\int }}\underset{\mathrm{0}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}\underset{\mathrm{0}}{\overset{{}^{\mathrm{2}\mathit{\pi }}}{\int }}\underset{\mathrm{0}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}{\mathbf{R}}_{\mathrm{g}}({\mathit{\theta }}^{\prime },{\mathit{\phi }}^{\prime };\mathit{\theta },\mathit{\phi })\cos {\mathit{\theta }}^{\prime }\sin {\mathit{\theta }}^{\prime }\cos \mathit{\theta }\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\phi }\mathrm{d}{\mathit{\theta }}^{\prime }\mathrm{d}{\mathit{\phi }}^{\prime }.\end{array}$$

The blue-sky albedo is a linear combination of A_bs and A_ws weighted by the fraction of diffuse skylight. The use of either A_bs, A_ws or A_LER in different applications depends on the assumptions of each parameter and the particular application. In the NIR, the diffuse component in the radiation field is much smaller than the direct component so the use of the A_bs is justified because it assumes only direct light. A_ws might be more suitable for applications in the UV/visible spectral range, where the diffuse component of the incident light may be of comparable size to the direct component. In any case, MODIS visible and NIR A_ws and A_bs do not differ on average by more than 5 % in summer (Oleson et al., 2003). A_ws is constant with solar incident direction, so its use is valid as a A_LER but it accounts for some surface BRDF effects (Eq. 3).

Figure 3 Sketch of the surface to top-of-atmosphere system with zenith and azimuth angles that define incident (θ^′, φ^′) and reflected (θ, φ ) directions of light. Direct incident solar light is described by (θ₀, φ₀). A_LER is the Lambertian surface albedo and R_g the surface BRDF.

Several models have been developed to describe surface BRDF (Wanner et al., 1995). These are either physical, empirical, or semi-empirical models. Physical models are constructed using the laws of physics to explicitly describe the processes that lead to the anisotropic behaviour of surface reflectance. Empirical models characterise the BRDF using mathematical functions that are suitable to describe the observed surface reflectance. Semi-empirical models describe the surface BRDF as a weighted sum of empirical functions derived from physical approximations.

Semi-empirical models are commonly used for global surface BRDF characterisation using remote sensing instruments. In these models, surface reflectance is represented as a linear combination of different terms (the kernels) that characterise different types of scattering that lead to the directional signatures on the reflectance (Roujean et al., 1992). Typically these terms consist of an isotropic term, a volume scattering term, and a geometric scattering term. The weights of the kernels cannot be directly interpreted as physical characteristics from the reflecting surface but just as a first-order approximation of the structure of the surface BRDF (Gao et al., 2003).

In the semi-empirical BRDF Ross Thick – Li Sparse (hereinafter Ross–Li) kernel-driven model, the surface reflectance is expressed as the sum of an isotropic term and two kernels (K_i) that depend on incident zenith angle (θ^′), viewing zenith angle (θ), and relative azimuth angle ($\mathit{\phi }-{\mathit{\phi }}^{\prime }$):

$$\begin{array}{ll}\text{(4)}& {\displaystyle }{\mathbf{R}}_{\mathrm{g}}(\mathit{\theta },{\mathit{\theta }}^{\prime },\mathit{\phi }-{\mathit{\phi }}^{\prime },\mathit{\lambda })=& {\displaystyle }{f}_{\mathrm{iso}}\left(\mathit{\lambda }\right)+f_{\mathrm{vol}}\left(\mathit{\lambda }\right)K_{\mathrm{vol}}(\mathit{\theta },{\mathit{\theta }}^{\prime },\mathit{\phi }-{\mathit{\phi }}^{\prime }){\displaystyle }& {\displaystyle }+f_{\mathrm{geo}}\left(\mathit{\lambda }\right)K_{\mathrm{geo}}(\mathit{\theta },{\mathit{\theta }}^{\prime },\mathit{\phi }-{\mathit{\phi }}^{\prime }).\end{array}$$

In Eq. (4) K_i are the kernels from the semi-empirical Ross–Li model that describe the three basic scattering types and f_i are the surface BRDF parameters that are retrieved from surface reflectance observations from satellite measurements (e.g. MODIS). The surface reflectance cloud-free observations used to obtain the surface BRDF parameters are corrected for absorption and scattering by atmospheric gases, aerosols, and thin clouds (Vermote et al., 1997). Improvements in the atmospheric correction scheme include as much information as possible derived from the satellite itself (e.g. MODIS aerosol optical thickness) (Vermote and Kotchenova, 2008).

The isotropy parameter (f_iso) represents the isotropic scattering from a nadir incident and nadir view position. The volumetric scattering is represented by the Ross-Thick kernel, K_vol. It is derived in a single-scattering approximation from radiative transfer theory for a thick homogeneous layer of small scatterers, with equal reflectance and transmittance (Roujean et al., 1992). To account for the reflectance peak in the backward-scattering direction (i.e. the hotspot effect), we include the modification on the volumetric kernel by Maignan et al. (2004). The geometric scattering is represented by the Li–Sparse kernel, K_geo. For this case, the scene is assumed to contain sparse objects that cast perfectly black shadows with sunlit and shaded portions of the ground and crown contributing to the modelled reflectance of the scene (Li and Strahler, 1986). The exact formulae for these kernels are summarised in Appendix A.

3.2 Surface BRDF implementation in DAK

The radiative transfer model DAK (Doubling–Adding KNMI, Lorente et al., 2017; Stammes et al., 1989; de Haan et al., 1987) is used in the GOME-2 and OMI cloud retrievals (FRESCO and OMCLDO2) and in the DOMINO NO₂ retrieval. Originally DAK only considered Lambertian surfaces. To account for the surface reflectance anisotropy in DAK, we have implemented the Ross–Li kernel-driven model. We chose this model for consistency with the retrieval algorithm of the MODIS BRDF/albedo product that we will use in our simulations. The MODIS satellite provides a reliable surface BRDF product and its resolution is suitable to capture surface anisotropy variations for OMI and GOME-2 resolution. MODIS BRDF/albedo products have been successfully used in different fields of atmospheric and climate science such as analysis of radiative forcing due to vegetation change (Myhre et al., 2005) and assessment of land surface albedo in global climate models (Oleson et al., 2003).

After the implementation of the Ross–Li surface BRDF model in DAK, the surface reflectance matrix R_g now contains the full reflection properties of the surface and substitutes the constant isotropic value used for the Lambertian case. This matrix is filled with the surface reflectance calculated with the Ross–Li BRDF model via Eq. (4) as a function of $\mathit{\theta },{\mathit{\theta }}^{\prime }$ for a specific $\mathit{\phi }-{\mathit{\phi }}^{\prime }$. For the Lambertian case, the matrix only contains the (1,1) element, which is the value of the surface albedo. We neglect polarisation in the BRDF. In the Doubling–Adding method for radiative transfer calculation, all the matrices (scattering, reflection, and transmission matrices) are expanded in a Fourier series for the integration over $\mathit{\phi }-{\mathit{\phi }}^{\prime }$ following the approach in de Haan et al. (1987). For each Fourier term the Doubling–Adding procedure is applied separately, including the addition of the surface. The mth Fourier coefficient matrix for the surface reflectance matrix is obtained from the relation:

$$\begin{array}{ll}\text{(5)}& {\displaystyle }& {\displaystyle }{\mathbf{R}}_{\mathrm{g}}^m(\mathit{\mu },{\mathit{\mu }}^{\prime })={\displaystyle }& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}d(\mathit{\phi }-{\mathit{\phi }}^{\prime })\cos m(\mathit{\phi }-{\mathit{\phi }}^{\prime }){\mathbf{R}}_{\mathrm{g}}(\mathit{\mu },{\mathit{\mu }}^{\prime },\mathit{\phi }-{\mathit{\phi }}^{\prime }),\end{array}$$

where μ, μ^′ are the cosines of the zenith angles in the scattered and incident direction ($\mathit{\theta },{\mathit{\theta }}^{\prime }$) respectively and $\mathit{\phi }-{\mathit{\phi }}^{\prime }$ is the difference between the scattered and incident azimuth angles.

Figure 4 Surface reflectance modelled with the Ross–Li surface BRDF model with parameters from MODIS band 3 (459–479 nm) representing a vegetated surface over Amazonia (a) as function of viewing zenith angle in the principal plane ($\mathit{\phi }-{\mathit{\phi }}^{\prime }=\mathrm{180}{}^{\circ }$ for negative θ and $\mathit{\phi }-{\mathit{\phi }}^{\prime }=\mathrm{0}{}^{\circ }$ for positive θ) and (b) in a polar plot with $\mathit{\phi }-{\mathit{\phi }}^{\prime }$ along the azimuth axis and θ along the polar axis for ${\mathit{\theta }}^{\prime }=\mathrm{30}{}^{\circ }$.

The coefficients of the Fourier expansion (Eq. 5) are calculated with the Gauss–Legendre quadrature integration method. It is possible to apply this method because the surface reflection matrix R_g is known at a certain number of division points $\mathit{\phi }-{\mathit{\phi }}^{\prime }$. The number of Fourier terms (m) and Gaussian points for azimuth ($\mathit{\phi }-{\mathit{\phi }}^{\prime }$) integration needed to resolve the surface BRDF shape depends on the illumination and viewing geometry. For geometries close to the hotspot region in the backward-scattering direction, the number of Fourier terms and Gaussian points needed to reproduce the original BRDF increases significantly with respect to geometries outside the hotspot. In order to reach an accuracy of 10⁻³ (difference between the original BRDF and the reconstructed BRDF with the Fourier expansion) over the hotspot region, 720 Gaussian points are needed for the azimuth integration, and 300 Fourier terms. Outside the hotspot region, using 60 Gaussian points and 30 Fourier terms in DAK reproduces the original surface reflectance values with an accuracy higher than 10⁻⁵. In the final implementation of the surface BRDF in the DAK RTM and in order to have an optimal simulation time, we used 100 Fourier terms and 360 Gaussian points for $\mathit{\phi }-{\mathit{\phi }}^{\prime }$, also over the hotspot. The overall accuracy obtained with these numbers is within the errors of the radiative transfer modelling for application in satellite retrievals (Lorente et al., 2017).

One of the disadvantages of empirical models is that they depend on observations to derive the parameters f_i that describe the surface BRDF. Kernel-based semi-empirical models like the Ross–Li model implemented in DAK only describe the surface reflection accurately for the range of illumination and viewing geometries of the measurements from which they have been derived (Litvinov et al., 2011). The geometries for which the semi-empirical models are valid are thus limited by extreme geometries of instruments like MODIS, which are typically 60–70^∘ for θ and θ^′. This limit overlaps well with the θ, θ^′ values in OMI and GOME-2 measurements for which clouds and trace gas columns are retrieved. For geometries outside the range of MODIS measurement geometries, surface reflectance variations need to be extrapolated.

In radiative transfer modelling with the Doubling–Adding method, in order to calculate the radiation field correctly, the reflectance and transmittance values are needed for all $\mathit{\theta },{\mathit{\theta }}^{\prime }$ between 0 and 90^∘ to be integrated over the complete hemisphere. However, after implementing the Ross–Li BRDF model in DAK, surface reflectance was negative or too high for some combinations of extreme geometries. In order to avoid these unphysical values, we tried various ways of extrapolating values from valid (MODIS-range) angles to the more extreme angles. The different methods did not affect TOA reflectance and albedo values by more than 1 %. Finally, we constrained the surface reflectance to the range [0,1] as in Eq. (6). This is reasonable as negative surface reflectance values are not physically valid, and in nature, surfaces with reflectance higher than 1 do not usually occur (except for surfaces covered by snow or ice). Therefore,

$$\begin{array}{}\text{(6)}& \begin{array}{lll}\mathrm{If}{\mathbf{R}}_{\mathrm{g}}({\mathit{\mu }}_{\mathrm{i}},{\mathit{\mu }}_{\mathrm{i}}^{\prime },{\mathit{\phi }}_{\mathrm{i}}-{\mathit{\phi }}_{\mathrm{i}}^{\prime })<\mathrm{0}& \mathrm{then}& {\mathbf{R}}_{\mathrm{g}}({\mathit{\mu }}_{\mathrm{i}},{\mathit{\mu }}_{\mathrm{i}}^{\prime },{\mathit{\phi }}_{\mathrm{i}}-{\mathit{\phi }}_{\mathrm{i}}^{\prime })=\mathrm{0}\\ \mathrm{If}{\mathbf{R}}_{\mathrm{g}}({\mathit{\mu }}_{\mathrm{i}},{\mathit{\mu }}_{\mathrm{i}}^{\prime },{\mathit{\phi }}_{\mathrm{i}}-{\mathit{\phi }}_{\mathrm{i}}^{\prime })>\mathrm{1}& \mathrm{then}& {\mathbf{R}}_{\mathrm{g}}({\mathit{\mu }}_{\mathrm{i}},{\mathit{\mu }}_{\mathrm{i}}^{\prime },{\mathit{\phi }}_{\mathrm{i}}-{\mathit{\phi }}_{\mathrm{i}}^{\prime })=\mathrm{1}.\end{array}\end{array}$$

Figure 4 is an example of the surface reflectance computed with the Ross–Li BRDF model for a surface with parameters (f_iso, f_vol, $f_{\mathrm{geo}})=(\mathrm{0.0399},\mathrm{0.0245},\mathrm{0.0072})$. This combination of f_i values are the spatially averaged parameters from MODIS (BRDF/albedo product MCD43A1) band 3 (459–479 nm) over Amazonia (lat: 5^∘ N–10^∘ S, long: 60–70^∘ W) for March 2008. This representation of a vegetated surface will be used in all the simulations in Sect. 3. The backward-scattering direction corresponds to values of $\mathit{\phi }-{\mathit{\phi }}^{\prime }=\mathrm{180}{}^{\circ }$ in the polar plot (Fig. 4b) and negative values of θ in the principal plane plot (Fig. 4a). In the backward-scattering direction, the surface reflectance is 2 times higher than in the forward-scattering direction ($\mathit{\phi }-{\mathit{\phi }}^{\prime }=\mathrm{0}{}^{\circ }$ in Fig. 4b and positive θ in Fig. 4a). In the backward-scattering direction the hotspot is clearly visible when $\mathit{\theta }=(-)\mathrm{30}{}^{\circ }={\mathit{\theta }}^{\prime }$.

The surface reflectance dependence on geometry as shown in Fig. 4 does not exist when a constant isotropic albedo or surface LER is used. Using a constant albedo for this surface of 0.03 (the equivalent A_ws is 0.034) means that in the backward-scattering direction the surface reflectance is underestimated. On the contrary, in the forward direction the surface reflectance is overestimated if an isotropic surface albedo is used. This surface reflectance difference between forward–backward-scattering direction, thus qualitatively explains the east–west bias in the surface LER in Fig. 1. West measurements in GOME-2 correspond to the backward-scattering direction, and for this direction the surface reflectance is higher than in the forward-scattering direction, i.e. east measurement in GOME-2. How this affects effective cloud fractions and NO₂ AMFs from GOME-2 and OMI will be analysed in Sect. 4.

Figure 5 TOA reflectance at (a) 469 nm and (b) 645 nm simulated by DAK (blue), SCIATRAN (red), and LIDORT (green) with the Ross–Li surface BRDF model with parameters representing a vegetated surface over Amazonia as in Fig. 4 as a function of the viewing zenith angle in the principal plane ($\mathit{\phi }-{\mathit{\phi }}_{\mathrm{0}}=\mathrm{180}{}^{\circ }$ for negative θ and $\mathit{\phi }-{\mathit{\phi }}_{\mathrm{0}}=\mathrm{0}{}^{\circ }$ for positive θ), for ${\mathit{\theta }}_{\mathrm{0}}=\mathrm{30}{}^{\circ }$.

3.3 Evaluation of surface BRDF effects in DAK TOA reflectances

To evaluate the surface BRDF implementation in DAK we compare TOA reflectances with those simulated by other state-of-the-art radiative transfer models that include a description of surface BRDF effects. Both SCIATRAN (Rozanov et al., 2014) and LIDORT (Spurr, 2004) model the surface reflectance using the Ross–Li model. To minimise the differences due to factors other than the surface BRDF implementations itself, the settings of the three RTMs were as similar as possible. These settings include no polarisation, a plane-parallel standard midlatitude atmosphere and absorption by O₃, O₂–O₂, and NO₂.

We select two combinations of the Ross–Li BRDF parameters to model surface reflectance of different surfaces. We simulate TOA reflectances at two different wavelengths: 469 and 645 nm. These wavelengths correspond to the middle of the MODIS band 3 (459–479 nm) and band 1 (620–670 nm) respectively.

Figure 5 shows the simulated TOA reflectance by DAK, SCIATRAN and LIDORT as a function of VZA along the principal plane at 469 and 645 nm. The agreement between the models is within 0.5 % for geometries outside the hotspot region, and DAK and SCIATRAN agree within 1 % over the hotspot. With this simple validation we ensure that our surface BRDF implementation in DAK is correct. The viewing geometry dependency of the TOA reflectance is the effect of Rayleigh scattering of a clean atmosphere. By comparing Fig. 4a with Fig. 5 we see that the effect of surface reflectance anisotropy is dampened by scattering in the atmosphere but the hotspot effect is still visible at both wavelengths after the radiation has passed through the atmosphere. Figure 5 also shows that the hotspot effect is less prominent in TOA reflectance at 469 nm compared to 645 nm, where the atmosphere is more transparent. At 469 nm Rayleigh scattering is stronger than at 645 nm, reducing the effects of surface anisotropy on the TOA reflectance.

We now compare TOA reflectance simulated with surface BRDF and TOA reflectance simulated assuming a Lambertian surface at 469 nm (Fig. 6a, b) and 758 nm (Fig. 6c, d). For simulations at 758 nm we use surface BRDF parameters from MODIS band 2 (841–876 nm) to account for the increase in surface reflectivity near 700 nm (the red edge, e.g. Tilstra et al., 2017). To test the representativeness of band 2 at 758 nm, we scaled the parameters from band 3 (459–479 nm) using the ratio of reflectance at 772 and 469 nm, and the differences with the parameters from MODIS band 2 were negligible. For the Lambertian case, we use the equivalent A_ws calculated using the surface BRDF model with Eq. (3). Figure 6b shows that, at 469 nm, the highest absolute differences are around the hotspot region in the backward-scattering direction, where TOA reflectance with BRDF is up to 15 % higher than the TOA reflectance for the Lambertian case. In the forward-scattering direction, TOA reflectance with surface BRDF is 7 % lower than for the Lambertian surface. At 758 nm the effect of surface reflectance anisotropy is much stronger than at 469 nm (Fig. 6c). The highest absolute differences at 758 nm (Fig. 6d) are in the hotspot region in the backward-scattering direction (up to 30 %) and for very extreme angles ($\mathit{\theta }>\mathrm{85}{}^{\circ }$) in the forward-scattering direction.

These results are consistent with those shown in Figs. 1 and 2: for GOME-2A west measurements (backward-scattering direction), TOA reflectances with surface BRDF are higher than for the Lambertian surface. If these differences are not accounted for in the cloud retrieval, cloud fractions will be biased high in the west measurements to match the measured TOA reflectance. For east measurements in the forward-scattering direction, TOA reflectance with surface BRDF is lower than for the Lambertian surface. Cloud fractions retrieved with surface LER will be biased low to match the measured TOA reflectance. Results from Fig. 6 underline that surface BRDF effects in retrieved cloud fractions are stronger for FRESCO at 758 nm than for OMCLDO2 at 477 nm. Our results also show that the error in TOA reflectances due to the use of a Lambertian albedo is substantial, but its magnitude is highly dependent on the spectral and geometrical characteristics of the measurements: effects are stronger at 758 nm and around the hotspot region in the backward-scattering direction.

Figure 6 (a, c) TOA reflectance as a function of viewing zenith angle simulated by DAK at 469 and 758 nm with a Lambertian surface (blue line) and with surface BRDF (green line) in the principal plane ($\mathit{\phi }-{\mathit{\phi }}_{\mathrm{0}}=\mathrm{180}{}^{\circ }$ for negative θ and $\mathit{\phi }-{\mathit{\phi }}_{\mathrm{0}}=\mathrm{0}{}^{\circ }$ for positive θ). (b, d) Absolute differences between TOA reflectance with surface BRDF and with a Lambertian surface at 469 and 758 nm. Surface BRDF parameters represent a vegetated surface over Amazonia at 469 nm (f_iso, f_vol, $f_{\mathrm{geo}})=(\mathrm{0.0399},\mathrm{0.0245},\mathrm{0.0072})$ and 758 nm (f_iso, f_vol, $f_{\mathrm{geo}})=(\mathrm{0.36},\mathrm{0.24},\mathrm{0.03})$. Note the different scales.

4 Role of surface BRDF in cloud retrievals

4.1 Synthetic cloud fraction retrieval

In Sect. 2 we showed that there is an east–west bias in the retrieved cloud fractions from GOME-2 and OMI. Effective cloud fractions in FRESCO and OMCLDO2 are retrieved as follows (Stammes et al., 2008; Veefkind et al., 2016):

$$\begin{array}{}\text{(7)}& c_{\mathrm{eff}}={\displaystyle \frac{R_{\mathrm{meas}}-R_{\mathrm{cr}}}{R_{\mathrm{cd}}-R_{\mathrm{cr}}}}.\end{array}$$

R_meas is the TOA reflectance measured by the satellite instrument, R_cr is the simulated clear-sky TOA reflectance, and R_cd is the simulated cloudy-sky TOA reflectance assuming that the cloud is a Lambertian reflector with a fixed albedo of 0.8. In the current versions of FRESCO (v7) and OMCLDO2 (v2.0) cloud retrievals, the simulated clear-sky TOA reflectances in Eq. (7) assume that the surface is Lambertian. Due to this assumption, any surface anisotropy signal in the measured TOA reflectance (R_meas) is neglected and, consequently, ends up in the retrieved effective cloud fraction.

Table 1 Settings for Lambertian and BRDF c_eff simulations in Sect. 4.1. Inverse model for Lambertian c_eff reproduces current O₂–O₂ and FRESCO+ retrievals, with Lambertian surface and Lambertian cloud. Inverse model for BRDF c_eff reproduces the retrieval accounting for surface BRDF effects.

Figure 7 Simulated effective cloud fraction at (a) 477 nm and (b) 758 nm as a function of viewing zenith angle for different geometric cloud fractions along the principal plane ($\mathit{\phi }-{\mathit{\phi }}_{\mathrm{0}}=\mathrm{180}{}^{\circ }$ for negative θ and $\mathit{\phi }-{\mathit{\phi }}_{\mathrm{0}}=\mathrm{0}{}^{\circ }$ for positive θ). ${\mathit{\theta }}_{\mathrm{0}}=\mathrm{30}{}^{\circ }$ and c_geo=0, 0.05, 0.2, and 0.5. At 477 nm: BRDF parameters (solid line) (f_iso, f_vol, $f_{\mathrm{geo}})=(\mathrm{0.03},\mathrm{0.02},\mathrm{0.01})$ and for Lambertian surface A_ws=0.0217 (dashed line). At 758 nm: BRDF parameters (solid line) (f_iso, f_vol, $f_{\mathrm{geo}})=(\mathrm{0.4},\mathrm{0.25},\mathrm{0.08})$ and for Lambertian surface A_ws=0.337 (dashed line).

To improve our understanding of how surface reflectance anisotropy influences the retrieval of cloud fractions, we use the forward model DAK to approximate R_meas by simulating the TOA reflectance for a scene with a Henyey–Greenstein cloud and surface reflectance anisotropy. This resembles what the satellite would measure in a realistic cloudy scene. We express R_meas as the sum of TOA reflectance of the cloudy and the clear parts of the scene, weighted by a geometric cloud fraction c_geo (independent pixel approximation):

$$\begin{array}{}\text{(8)}& R_{\mathrm{meas}}\approx R_{\mathrm{sim}}=c_{\mathrm{geo}}\cdot R_{\mathrm{cd}}+(\mathrm{1}-c_{\mathrm{geo}})R_{\mathrm{cr}}.\end{array}$$

The effective cloud fraction is the part of the pixel that the Lambertian cloud has to occupy to match the observed reflectance. The geometric cloud fraction is the part of the pixel that is covered by the “true” cloud (Stammes et al., 2008).

The settings of the simulations are summarised in Table 1. We simulate TOA reflectances and cloud fractions in the spectral regions where cloud fractions are calculated in the cloud algorithms. The Ross–Li BRDF parameters are from a climatology created by the QA4ECV Land Group at the Mullard Space Science Laboratory (University College London). This data set consists of daily BRDF parameters collected from 16 years of MODIS measurements (Strahler et al., 1999; MCD43A1, 2015) from 2000 to 2016 (QA4ECV-WP4, 2016). Parameters from band 3 (459–479 nm) are representative for simulations in the O₂–O₂ absorption band and parameters from band 2 (841–876 nm) for simulations in the O₂–A band (see Sect. 3.2). Monthly averaged parameters from this data set over Amazonia are shown in Figs. S2, S3 in the Supplement.

We calculate cloud fractions using Eq. (7) with R_cr simulated with a Lambertian surface consistent with the current O₂–O₂ and FRESCO+ retrievals (hereafter Lambertian c_eff) and by accounting for surface BRDF effects (hereafter BRDF c_eff). The Lambertian cloud is located at the same pressure level as the Henyey–Greenstein cloud so we can isolate surface BRDF effects on cloud fraction only (see settings in Table 1).

Figure 7 shows that cloud fractions accounting for surface BRDF effects (dashed lines) depend only weakly on geometry, whereas cloud fractions with a Lambertian surface (solid lines) are higher in the backward-scattering direction, at both 477 and 758 nm and especially for the lowest c_geo. At 758 nm this is true even for a geometric cloud fraction as high as 0.5. In the backward-scattering direction, surface reflectance is higher than reflectance by a Lambertian surface. Therefore clear-sky TOA reflectance with a Lambertian surface cannot explain the higher-simulated reflectance in Eq. (7), which results in high Lambertian c_eff.

Surface BRDF effects are more important for small cloud fractions, and less importance for large cloud fractions. For cloudy pixels the effect of surface reflectance anisotropy vanishes because the scattering by the cloud dominates in the measured reflectance. Figure 7 shows that for large c_geo, the viewing zenith angle dependency of the effective cloud fractions is due to the Henyey–Greenstein scattering cloud, which gives relatively higher scattering in the forward direction.

At 477 nm, surface BRDF effects on cloud fractions are less evident than at 758 nm. In the visible spectral region, surface BRDF effects are suppressed by Rayleigh scattering smoothing out the surface anisotropy effects on TOA reflectances. For lower geometric cloud fraction, Lambertian cloud fractions are moderately higher (by 0.05) in the backward-scattering direction than in the forward-scattering direction. These findings underscore the relevance of accounting for surface BRDF effects because measurements with small cloud fractions are most sensitive to pollution in the lower troposphere.

Differences in Lambertian cloud fractions between backward- and forward-scattering directions at 758 nm are on average 0.35. At 477 nm, the differences amount to 0.1, depending on the surface and the geometry. This is consistent with the observed bias in FRESCO and OMCLDO2 cloud fractions shown in Fig. 2. The absence of a backward–forward-scattering dependency in the BRDF c_eff implies that accounting for surface BRDF effects will reduce the east–west across-track bias in retrieved cloud fractions.

4.2 GOME-2A cloud fraction simulations

We simulate Lambertian and BRDF cloud fractions for GOME-2A measurements over Amazonia in March 2008. We use the exact illumination and viewing geometry (θ,θ₀, φ−φ₀) of each individual measurement and co-locate MODIS pixels with the GOME-2A pixel centre to obtain the surface BRDF parameters over the scene. For Lambertian cloud fractions, we use the GOME-2A surface LER value from Tilstra et al. (2017). To simulate the measured reflectance, we assume a geometric cloud fraction distribution with an area-wide average of $\stackrel{\mathrm{‾}}{c_{\mathrm{geo}}}=\mathrm{0.33}$. Figure 8a, d show the $\stackrel{\mathrm{‾}}{c_{\mathrm{geo}}}$ distribution for east and west measurements respectively.

Figure 8 shows simulated Lambertian (panels b, e) and BRDF (panels c, f) effective cloud fractions for east and west GOME-2A measurements. For east measurements, both Lambertian ($\stackrel{\mathrm{‾}}{c_{\mathrm{eff}}}=\mathrm{0.32}$) and BRDF ($\stackrel{\mathrm{‾}}{c_{\mathrm{eff}}}=\mathrm{0.28}$) cloud fraction simulations capture the original geometric cloud fraction distribution and absolute values. This means that surface BRDF effects are weaker in the east, consistent with the smaller surface LER climatology bias for measurements on the east part of the orbit (Fig. 1). For west measurements, the original distribution is reproduced well in the BRDF cloud fractions ($\stackrel{\mathrm{‾}}{c_{\mathrm{eff}}}=\mathrm{0.29}$), but the Lambertian values show a significant overestimation ($\stackrel{\mathrm{‾}}{c_{\mathrm{eff}}}=\mathrm{0.51}$). A box plot of these distributions is shown in Fig. S4 in the Supplement. The east–west bias in the Lambertian simulation (0.2) is considerably reduced (0.01) for the BRDF cloud fractions: the across-track bias in the FRESCO data is a direct consequence of neglecting surface BRDF effects. We conclude that accounting for these surface BRDF effects can largely solve the bias in cloud fractions measured in the backward-scattering regime over Amazonia. Although we have not made an analysis of the surface BRDF effects in a complete retrieval, the biases in cloud fraction found over other regions will probably be reduced after accounting for surface BRDF effects.

Figure 8 Geometric cloud fraction distributions for (a) east and (d) west GOME-2A measurements. (b, e) Lambertian and (c, f) BRDF cloud fraction simulations for east (a–c) and west (d–f) GOME-2A measurements. Plots show averaged cloud fractions in a $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ grid over Amazonia (lat: 5^∘ N–10^∘ S, long: 60–70^∘ W) for March 2008.

5 Role of surface BRDF in NO₂ retrievals

Here we investigate the effects of accounting for surface reflectance anisotropy on tropospheric NO₂ column retrievals under clear-sky and partly cloudy conditions. We calculate tropospheric AMFs with surface BRDF and with a Lambertian surface, including the surface BRDF effects on retrieved cloud fractions (Sect. 4) for GOME-2A measurements over Amazonia and over France.

In satellite retrievals of trace gases, an air mass factor (M) is used to convert the slant column density (N_s, SCD) from the measured reflectance spectra into the vertical column density (N_v, VCD):

$$\begin{array}{}\text{(9)}& M(R_{\mathrm{s}},P_{\mathrm{s}},f_{\mathrm{cd}},P_{\mathrm{cd}},{\mathit{x}}_{\mathrm{a}},\mathit{\theta },{\mathit{\theta }}_{\mathrm{0}},\mathit{\phi }-{\mathit{\phi }}_{\mathrm{0}})={\displaystyle \frac{N_{\mathrm{s}}}{N_{\mathrm{v}}}}.\end{array}$$

The AMF depends on the surface reflectance (R_s), surface pressure (P_s), cloud fraction, and pressure (f_cd,P_cd), a priori NO₂ profile (x_a) and measurement geometry ($\mathit{\theta },{\mathit{\theta }}_{\mathrm{0}},\mathit{\phi }-{\mathit{\phi }}_{\mathrm{0}}$).

Here, tropospheric NO₂ AMFs are calculated by differencing the logarithm of simulated TOA reflectances with and without trace gas in the troposphere divided by the absorption optical thickness of the gas τ_gas:

$$\begin{array}{}\text{(10)}& M=-{\displaystyle \frac{\ln R\left({\mathit{\tau }}_{\mathrm{gas}}\right)-\ln R({\mathit{\tau }}_{\mathrm{gas}}=\mathrm{0})}{{\mathit{\tau }}_{\mathrm{gas}}}}.\end{array}$$

Table 2 Settings for the Lambertian and BRDF tropospheric NO₂ AMF calculations shown in Fig. 9.

AMF is directly affected by the assumption of a Lambertian surface instead of an anisotropic surface in the simulated TOA reflectance. In addition, AMFs are indirectly affected by the cloud radiance fraction (f_cd) used to correct for residual clouds, in which calculation a Lambertian surface is assumed as well. To account for the presence of clouds, we use the independent pixel approximation (IPA; see also Eq. 8), which consists of calculating the total AMF for a partly cloudy scene as a linear combination of cloudy (M_cd) and clear (M_cr) components of the AMF, weighted by the cloud radiance fraction w:

$$\begin{array}{}\text{(11)}& M=w{M}_{\mathrm{cd}}+(\mathrm{1}-w)M_{\mathrm{cr}}.\end{array}$$

The use of a Lambertian surface thus influences the AMF directly via M_cr and indirectly via w:

$$\begin{array}{}\text{(12)}& w={\displaystyle \frac{c_{\mathrm{eff}}{R}_{\mathrm{cd}}}{c_{\mathrm{eff}}{R}_{\mathrm{cd}}+(\mathrm{1}-c_{\mathrm{eff}})R_{\mathrm{cr}}}},\end{array}$$

where c_eff is the cloud fraction, and R_cd and R_cr are the radiances for totally cloudy and clear-sky scenes. Cloud radiance fraction depends on the Lambertian surface assumption via R_cr and c_eff.

5.1 BRDF effects on tropospheric NO₂ air mass factors

We calculate AMFs with Eq. (10) by simulating TOA reflectances with surface BRDF and with a Lambertian surface. The settings of the simulations are summarised in Table 2. Based on our analysis (Sect. 4.1), we include a change of ±0.05 over the Lambertian cloud fraction with a decrease of −0.05 in the backward-scattering direction and an increase of +0.05 in the forward-scattering direction to quantify how the surface BRDF effects on clouds propagate to the final AMF.

Figure 9a shows surface BRDF effects on M_cr for a moderately polluted troposphere as a function of VZA along the principal plane. In the backward-scattering direction, BRDF M_cr is higher by 5–20 %. The higher surface BRDF reflectance and TOA reflectance makes the retrieval more sensitive to the NO₂ in the boundary layer. In the forward-scattering direction, BRDF M_cr is lower by 5–15 %.

Figure 9b shows surface BRDF effects on w and Fig. 9c shows the combined effect on total tropospheric M of changes in M_cr and in w in a partly cloudy scene with a c_eff(Lamb.)=0.1. In the backward-scattering direction, tropospheric M is 9–30 % higher when accounting for surface BRDF effects. The decrease in c_eff (and hence w) makes the retrieval more sensitive to the NO₂ below the cloud. In the forward-scattering direction, tropospheric M is 14 %–22 % lower when accounting for surface BRDF effects because of the stronger screening effect by the higher BRDF c_eff. The surface BRDF effect on w enhances the effect on the clear-sky AMF by up to 10 % in both the forward and backward-scattering directions.

Figure 9 (a) Clear-sky tropospheric NO₂ AMF, (b) cloud radiance fraction, and (c) total tropospheric NO₂ AMF computed with surface BRDF (green) and with a Lambertian surface (blue) and their relative difference (right axis) as a function of viewing zenith angle in the principal plane ($\mathit{\phi }-{\mathit{\phi }}_{\mathrm{0}}=\mathrm{180}{}^{\circ }$ for negative θ and $\mathit{\phi }-{\mathit{\phi }}_{\mathrm{0}}=\mathrm{0}{}^{\circ }$ for positive θ), for ${\mathit{\theta }}_{\mathrm{0}}=\mathrm{30}{}^{\circ }$. P_cd=850 hPa, Lambertian c_eff=0.1 and BRDF $c_{\mathrm{eff}}=\mathrm{0.10}\pm \mathrm{0.05}$. Surface BRDF parameters are (f_iso, f_vol, $f_{\mathrm{geo}})=(\mathrm{0.04},\mathrm{0.03},\mathrm{0.008})$ and A_ws=0.036 for the Lambertian surface. Troposphere is moderately polluted ($N_{\mathrm{v},\mathrm{trop}}=\mathrm{4}\times {\mathrm{10}}^{\mathrm{15}}$ molec cm⁻²).

Figure 10 shows surface BRDF effects on total tropospheric M in partly cloudy scenes with increasing cloudiness for the specific combination of (θ, ${\mathit{\theta }}_{\mathrm{0}})=(\mathrm{45}{}^{\circ }$, 30^∘) in the principal plane, representative of typical GOME-2 or OMI measurements. For a polluted troposphere (represented by stars), the total effect of accounting for surface BRDF is to increase M in the backward-scattering direction (Fig. 10a) and to reduce M in the forward-scattering direction (Fig. 10b). Average relative differences of about 15 % for M_cr increase up to 25 %–40 % for low cloud fractions (below 0.1) (see Fig. S5 in the Supplement). The effect decreases for higher cloud fractions, with relative differences below 10 % for cloud fractions higher than 0.5. For a lower bias in c_eff (e.g. 0.02), the effect on M reduces to about 15 %–20 %, which is still a considerable effect.

For an unpolluted troposphere, surface BRDF effects on M_cr are in the same direction as in the polluted case but smaller (−7 % to +11 %). In a partly cloudy scene, Lambertian and BRDF tropospheric M are similar within −3 % to 7 %. Because there is very little NO₂ below the cloud, the effect of a change in w counteracts (and largely cancels) the effect on the clear-sky AMF (circles in Fig. 10).

Figure 10 Total tropospheric NO₂ AMF as a function of cloud fraction in the (a) backward-scattering direction and (b) forward-scattering direction computed with surface BRDF (green) and a Lambertian surface (blue), for P_cd=850 hPa, (θ, ${\mathit{\theta }}_{\mathrm{0}})=(\mathrm{45}{}^{\circ }$, 30^∘), for a moderately polluted (stars) and unpolluted (circles) troposphere. BRDF parameters are (f_iso, f_vol, $f_{\mathrm{geo}})=(\mathrm{0.04},\mathrm{0.03},\mathrm{0.008})$ and A_ws=0.036 for the Lambertian surface.

Figure 11 shows surface BRDF effects on total tropospheric AMF as a function of cloud pressure, for a cloud fraction of 0.1. For cloud pressures higher than the 850 hPa assumed in Fig. 10, the contribution of surface BRDF effects to the change in M from the change in cloud fractions is dampened. There is a range of cloud pressures (in Fig. 11 between 900 and 950 hPa) for which the effects on M_cr and on cloud fraction compensate each other. For an even higher cloud pressure (e.g. 978 hPa), M_cd is larger than M_cr and the sign of the effect changes. In the backward-scattering we have lower BRDF AMFs and in the forward-scattering higher BRDF AMFs. In the unpolluted situations the differences are larger for higher cloud pressures.

Although this study does not address surface BRDF effects on cloud pressure, we did a preliminary analysis, applying a directional surface LER derived from GOME-2 in FRESCO+. The analysis shows that accounting for surface reflectance anisotropy effects reduces the cloud pressure by 40 hPa on average (with differences up to 120 hPa). This high bias in retrieved cloud pressure implies that the results shown for 850 hPa might be representative of the surface BRDF effects on AMFs for clouds currently retrieved at higher (biased) pressures.

Figure 11 Total tropospheric NO₂ AMF as a function of cloud fraction in the (a) backward-scattering direction and (b) forward-scattering direction computed with surface BRDF (green) and a Lambertian surface (blue), for P_cd=850 hPa, (θ, ${\mathit{\theta }}_{\mathrm{0}})=(\mathrm{45}{}^{\circ }$, 30^∘), for a moderately polluted (stars) and unpolluted (circles) troposphere. BRDF parameters are (f_iso, f_vol, $f_{\mathrm{geo}})=(\mathrm{0.04},\mathrm{0.03},\mathrm{0.008})$ and A_ws=0.036 for the Lambertian surface.

5.2 GOME-2A tropospheric NO₂ air mass factors

We calculate Lambertian and BRDF NO₂ tropospheric AMFs for the exact illumination and viewing geometries of GOME-2A measurements over Amazonia and over France for March 2008. This results in approximately 1300 clear-sky pixels analysed over Amazonia and 700 over France. We assume a moderately polluted atmosphere in every scene and collocate MODIS pixels with the GOME-2A pixel centre to obtain the surface BRDF parameters. For the Lambertian simulations, we use the surface LER value of each pixel from the GOME-2 climatology. We apply Lambertian and BRDF c_eff distributions from Sect. 4.2 (as in Fig. 8). This way we account for the calculated surface BRDF effects in cloud fraction instead of the average change of 0.05 assumed in the sensitivity analysis in Sect. 5.1.

Figure 12a shows that for east measurements (i.e. forward scattering), M_cr decreases on average by 18 % over Amazonia. Over France the decrease is on average 8 % (not shown). Figure 12b shows that for west measurements (i.e. backward scattering), M_cr increases on average 5 % over Amazonia and 7 % over France, consistent with our findings in Fig. 9a. The differences of 15 %–23 % found in this analysis agree with the reported differences of 10 %–20 % in Noguchi et al. (2014) and Zhou et al. (2010) for clear-sky AMFs. The higher BRDF AMFs in the upper-right corner of our study area correspond to a savanna ecosystem that is brighter than the dense rainforest. Because of the higher resolution of the MODIS BRDF data set, this feature is well captured and leads to higher AMFs both in the east and west measurements (see Fig. S2 in the Supplement). MODIS albedo is able to capture spatial variability on scales that satellites with coarser pixels cannot (Russell et al., 2011).

Figure 12 Relative differences between BRDF and Lambertian AMFs over Amazonia for March 2008: (a) clear-sky tropospheric NO₂ AMFs for east GOME-2A measurements (forward scatter), (b) clear-sky tropospheric NO₂ AMFs for west GOME-2A measurements (backward scatter), (c) total tropospheric NO₂ AMF for east GOME-2A measurements (forward scatter) and (d) total tropospheric NO₂ AMF for west GOME-2A measurements (backward scatter). Total AMFs are only shown for c_eff<0.5.

Figure 12c shows that for east measurements (i.e. forward scattering), total tropospheric AMFs including surface BRDF effects on cloud fractions are on average 10 % lower over Amazonia. The decrease is on average 7 % over France (not shown). For west measurements, (i.e. backward scattering), total AMFs are on average 16 % higher over both Amazonia (Fig. 12d) and France, illustrating that the M_cr effect is enhanced by the effect on the cloud fraction by 10 %–15 % on average. Vasilkov et al. (2017) also found increased differences in the tropospheric AMFs due to surface BRDF effects on cloud parameters, but as reported by Lin et al. (2014), there can also be compensating effects.

These results show that surface BRDF affects both clear-sky AMF and cloud radiance fractions, which in combination significantly affect total NO₂ AMFs. As shown in Figs. 10–12, the sign and magnitude of the surface BRDF effects show strong spatial variations and depend on cloud fraction and cloud pressure. In order to generalise the effects to a global retrieval, a full assessment including all possible retrieval conditions should be done. Over forested terrain, current tropospheric AMFs are likely underestimated in the backward-scattering regime and overestimated in the forward-scattering regime by up to 25 %–35 %, explained by systematic errors in M_cr and w. These results show that surface BRDF effects have to be included consistently in both cloud and trace gas retrievals.

6 Discussion and conclusions

We analysed the effects of surface reflectance anisotropy on the OMI and GOME-2 satellite retrievals of cloud parameters and tropospheric NO₂ columns that currently use Lambertian-equivalent reflectivity (LER) climatologies. These climatologies, and consequently retrieved cloud fractions, show substantial across-track biases over terrain with a strong BRDF directionality. Here we interpret these with the DAK radiative transfer model. A clear understanding of the reasons for the biases and how they propagate in the tropospheric NO₂ column retrieval is critical to improve cloud and trace gas retrieval algorithms for satellite sensors.

An important finding is that the LER climatologies slightly overestimate surface albedo for forward-scattering satellite-viewing geometries (eastern part of GOME-2 orbit) and highly underestimate the surface albedo for backward scattering viewing geometries (western part). The underestimation is as large as a factor of 2 over forested scenes in the near-infrared (772 nm). They are weaker but still relevant in the visible (494 nm), where surfaces are darker and Rayleigh scattering effects are stronger. This across-track bias in surface LER propagates into the cloud fraction retrievals: we find biases in cloud fractions of up to 50 % between backward-scattering and forward-scattering geometries in the GOME-2 FRESCO and 26 % in the OMI OMCLDO2 cloud algorithms. Time of day does not drive the importance of surface BRDF effects but specific viewing geometry and spectral range do.

To interpret the above biases, we extended the description of surface reflectance in DAK to include the geometrical surface reflecting properties via the bidirectional reflectance distribution function (BRDF) from the Ross–Li semi-empirical model. This allows DAK to simulate not only isotropic reflection at the surface but also the anisotropic contributions from volumetric (e.g. leaf scattering) and from geometric (e.g. shadow-casting) effects. We evaluated DAK top-of-atmosphere (TOA) reflectance simulations against other radiative transfer models and find agreement within 1 % between DAK and SCIATRAN, even within the hotspot backscatter reflectance peak. We then simulated TOA reflectances over vegetated scenes using BRDF parameters from a daily high-resolution database derived from 16-years of MODIS measurements recently developed within the QA4ECV project (QA4ECV-WP4, 2016). Our updated DAK simulations show considerably higher TOA reflectance levels for backward-scattering viewing geometries than those with isotropic surface reflection (LER) only. This strongly indicates that across-track biases in cloud fractions can be explained by the lack of a description of surface reflectance anisotropy in the FRESCO and OMCLDO2 algorithms.

Subsequent sensitivity tests indicated that accounting for surface reflectance anisotropy in the FRESCO and OMCLDO2 retrieval framework removes the bias in cloud fractions. A correct physical description of surface anisotropy is essential for FRESCO, because cloud properties are retrieved in the NIR spectral range (760–790 nm) where surface BRDF effects are stronger and the atmosphere is virtually transparent. It is also of high relevance for scenes with low cloud fractions, where trace gas retrievals are still sensitive to pollution close to the ground. A discussion on the validity of the Lambertian cloud model is beyond the scope of this study. Nevertheless, the cloud fraction dependency with VZA for cloudy scenes suggests that the use of a more realistic cloud model should be considered in future improvements of cloud retrievals.

The implications for NO₂ air mass factor (AMF) calculations are substantial. Total tropospheric NO₂ AMFs are calculated as the radiative cloud fraction-weighted sum of cloudy and clear-sky AMFs. For moderately polluted NO₂ and backward-scattering geometries, we find that clear-sky AMFs are up to 20 % higher and cloud radiance fractions up to 40 % lower if surface reflectance anisotropic effects are accounted for. The combined effect of these changes (with clouds located at 850 hPa) is that NO₂ AMFs in polluted situations increase by 25 %–30 % for backward-scattering geometries (and decrease by 25 %–35 % for forward-scattering geometries), stronger than the effect of either contribution alone.

An issue that was not addressed in this study is the role of aerosols. Noguchi et al. (2014) showed that scattering by aerosols generally dampens surface BRDF effects for clear-sky scenes. However, more research is needed to assess how specific aerosol characteristics (i.e. aerosol amount and type, vertical distribution relative to cloud) will affect cloud parameter retrievals and air mass factor calculations both in clear-sky and cloudy conditions.

We conclude that it is necessary to coherently account for surface reflectance anisotropy effects in retrievals of cloud properties and trace gases from UV/vis satellite sensors. Although this study does not apply surface BRDF to a complete global cloud and NO₂ retrieval, it shows that it has substantial effects on both cloud fractions and NO₂ AMFs. A number of recent studies have attempted to account for the effects of anisotropic reflectance on both cloud and NO₂ retrievals (Lin et al., 2014; Vasilkov et al., 2017), but a global assessment including the full range of possible retrieval conditions is still missing. An additional incentive to account for surface reflectance anisotropy is that the currently available LER climatologies (Kleipool et al., 2008; Tilstra et al., 2017) describe the spatial variation in albedo on a scale ($\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$–$\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$) coarser than the OMI or GOME-2 pixel itself (13×24 km²/80×40 km²). Using these coarse LER climatologies in AMF calculations degrades the intrinsic spatial resolution of the satellite retrievals, an issue that will be exacerbated for the recently launched TROPOMI instrument, with pixels as small as 3.5×7 km². A viable alternative to the current LER climatologies is provided by the MODIS-derived BRDF-parameters at a spatial resolution better than the GOME-2, OMI, and TROPOMI pixel sizes. MODIS Terra and Aqua are expected to last until 2025 and afterwards the Joint Polar Satellite System (JPSS) satellite constellation ensures continuity of land observations needed to produce surface BRDF data. Sentinel-3 could be employed to generate a BRDF similar to the one from the ESA GlobAlbedo broadband and the QA4ECV spectral albedo after some years of measurements. Another alternative is to make a directionally dependent LER database from TROPOMI once there is enough surface reflectance data acquired by the satellite itself.

Appendix A: Kernels for the Ross–Li BRDF model

We summarise here the expressions of the kernels implemented in DAK to model the surface reflectance anisotropy. We refer the reader to the original literature where these kernels were derived for more detailed information.

A1 Ross-Thick kernel

The expression of the Ross-Thick volumetric scattering kernel is (Roujean et al., 1992):

$$\begin{array}{}\text{(A1)}& K_{\mathrm{vol}}=\left[{\displaystyle \frac{(\mathit{\pi }/\mathrm{2}-\mathit{\xi })\cos \mathit{\xi }+\sin \mathit{\xi }}{\cos {\mathit{\theta }}^{\prime }+\cos \mathit{\theta }}}\right]\cdot \left[\mathrm{1}+(\mathrm{1}+\mathit{\xi }/{\mathit{\xi }}_{\mathrm{0}}{)}^{-\mathrm{1}}\right]-{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}.\end{array}$$

Here θ^′ and θ are the incident and reflected zenith angles respectively. ξ is the scattering angle defined as follows:

$$\begin{array}{}\text{(A2)}& \cos \mathit{\xi }=\cos \mathit{\theta }\cos {\mathit{\theta }}^{\prime }+\sin \mathit{\theta }\sin {\mathit{\theta }}^{\prime }\cos (\mathit{\phi }-{\mathit{\phi }}^{\prime }),\end{array}$$

where $\mathit{\phi }-{\mathit{\phi }}^{\prime }$ is the relative azimuth angle (reflected and incident azimuth difference).

The term in the second pair of squared bracket in Eq. (A1) is the modified part for the hotspot modelling, where ξ₀ is the hotspot characteristic angle (typically 1.5^∘). This characteristic angle can be related to the ratio of the size of the scattering element and the canopy vertical density (Maignan et al., 2004).

A2 Li–Sparse kernel

The expression of the Li–Sparse geometric scattering kernel is (Li and Strahler, 1986; Wanner et al., 1995):

$$\begin{array}{ll}\text{(A3)}& {\displaystyle }{K}_{\mathrm{geo}}=& {\displaystyle }O(\mathit{\theta },{\mathit{\theta }}^{\prime },\mathit{\phi }-{\mathit{\phi }}^{\prime })-\sec {\mathit{\theta }}^{*}-\sec {{\mathit{\theta }}^{\prime }}^{*}{\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{1}+\cos {\mathit{\xi }}^{*})\sec {\mathit{\theta }}^{*}\sec {{\mathit{\theta }}^{\prime }}^{*}.\end{array}$$

The different terms in the equation are as follows:

$$\begin{array}{}\text{(A4)}& O={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}(t-\sin t\cos t)(\sec {\mathit{\theta }}^{*}+\sec {{\mathit{\theta }}^{\prime }}^{*}),\end{array}$$

$$\begin{array}{ll}\text{(A5)}& {\displaystyle }& {\displaystyle }\cos t={\displaystyle \frac{h}{b}}{\displaystyle \frac{\sqrt{D^{\mathrm{2}}+\left(\tan {{\mathit{\theta }}^{\prime }}^{*}\tan {\mathit{\theta }}^{*}\sin \right(\mathit{\phi }-{\mathit{\phi }}^{\prime }){)}^{\mathrm{2}}}}{\sec {\mathit{\theta }}^{*}\sec {{\mathit{\theta }}^{\prime }}^{*}}},{\displaystyle }& {\displaystyle }D=\sqrt{{\tan }^{\mathrm{2}}{\mathit{\theta }}^{*}+{\tan }^{\mathrm{2}}{{\mathit{\theta }}^{\prime }}^{*}+\mathrm{2}\tan {\mathit{\theta }}^{*}\tan {{\mathit{\theta }}^{\prime }}^{*}\cos (\mathit{\phi }-{\mathit{\phi }}^{\prime })}.\end{array}$$

The angles with a star are equivalent angles to convert spheroid-like object to spheres:

$$\begin{array}{}\text{(A6)}& {\mathit{\theta }}^{*}={\tan }^{-\mathrm{1}}\left({\displaystyle \frac{b}{r}}\tan \mathit{\theta }\right),{{\mathit{\theta }}^{\prime }}^{*}={\tan }^{-\mathrm{1}}\left({\displaystyle \frac{b}{r}}\tan {\mathit{\theta }}^{\prime }\right).\end{array}$$

O (Eq. A4) is the overlap area between the shadow of illumination and the shadow of viewing projections on the ground. D is the distance between the centres of the scattering objects. The parameter t is used to parameterise the scattering objects spherically. This kernel it is not linear as it has two parameters $\frac{b}{r}$ and $\frac{h}{b}$ describing the shape and the relative height of the scattering objects, and in the simulations in this study are set to 1 and 2 respectively.