Directionally dependent Lambertian-equivalent reflectivity (DLER) of the Earth's surface measured by the GOME-2 satellite instruments

In this paper we introduce the new concept of directionally dependent Lambertian-equivalent reflectivity (DLER) of the Earth's surface retrieved from satellite observations. This surface DLER describes Lambertian (isotropic) surface reflection which is extended with a dependence on the satellite viewing geometry. We apply this concept to data of the GOME-2 satellite instruments to create a global database of the reflectivity of the Earth's surface, providing surface DLER for 26 wavelength bands between 328 and 772 nm as a function of the satellite viewing angle via a second-degree polynomial parameterisation. The resolution of the database grid is 0.25^∘ by 0.25^∘, but the real, intrinsic spatial resolution varies over the grid from 1.0^∘ by 1.0^∘ to 0.5^∘ by 0.5^∘ down to 0.25^∘ by 0.25^∘ by applying dynamic gridding techniques. The database is based on more than 10 years (2007–2018) of GOME-2 data from the MetOp-A and MetOp-B satellites.

The relation between DLER and bi-directional reflectance distribution function (BRDF) surface reflectance is studied using radiative transfer simulations. For the shorter wavelengths (λ<500 nm), there are significant differences between the two. For instance, at 463 nm the difference can go up to 6 % at 30^∘ solar zenith angle. The study also shows that, although DLER and BRDF surface reflectances have different properties, they are comparable for the longer wavelengths (λ>500 nm). Based on this outcome, the GOME-2 surface DLER is compared with MODIS surface BRDF data from MODIS band 1 (centred around 645 nm) using both case studies and global comparisons. The conclusion of this validation is that the GOME-2 DLER compares well to MODIS BRDF data and that it does so much better than the non-directional LER database. The DLER approach for describing surface reflectivity is therefore an important improvement over the standard isotropic (non-directional) LER approaches used in the past.

The GOME-2 surface DLER database can be used for the retrieval of atmospheric properties from GOME-2 and from previous satellite instruments like GOME and SCIAMACHY. It will also be used to support retrievals from the future Sentinel-5 UVNS (ultraviolet, visible, near-infrared, and short-wave infrared) satellite instrument.

1 Introduction

Most satellite retrievals of atmospheric composition require accurate information about the reflectivity of the Earth's surface to achieve accurate retrieval results. This includes the retrieval of trace gases, such as ozone (O₃), nitrogen dioxide (NO₂), bromine oxide (BrO), formaldehyde (CH₂O), water vapour (H₂O), carbon dioxide (CO₂), carbon monoxide (CO), and methane (CH₄), and of cloud and aerosol information. To date, many of these retrievals use Lambertian surface reflection in the radiative transfer calculations and, consequently, adopt the use of Lambertian (isotropic) surface albedo climatologies. Examples are the retrievals of NO₂ (e.g. Boersma et al., 2011; Bucsela et al., 2013), formaldehyde (CH₂O) (e.g. de Smedt et al., 2015; Hewson et al., 2015), and cloud products (e.g. Lelli et al., 2012; Veefkind et al., 2016). Although relying on Lambertian reflection is common practice, using a bi-directional reflectance distribution function (BRDF) (Nicodemus et al., 1992; Schaepman-Strub et al., 2006) to describe the surface reflectivity would be preferable. According to a recent study by Lorente et al. (2018), the simplification of using Lambertian surface reflection can lead to errors of a factor of 2 in the surface reflection for vegetated surfaces.

Recently, several different approaches have been introduced to address this issue. One example is the introduction of geometry-dependent surface Lambertian-equivalent reflectivity (GLER) (Vasilkov et al., 2017; Qin et al., 2019; Fasnacht et al., 2019). In the GLER approach, surface BRDF information from the MODIS surface BRDF database (Gao et al., 2005) is used to calculate Lambertian surface albedo at 466 nm for land-covered satellite footprints of the ozone monitoring instrument (OMI). For the footprints over water surfaces model calculations are used (Fasnacht et al., 2019). The result is a Lambertian surface albedo that is ready to be used in a radiative transfer code with Lambertian surface reflection, calculated for the exact scattering geometry of the OMI footprint and for the specific date of the OMI footprint. The advantage is that this Lambertian surface albedo is adjusted to the geometry of the observation, whereas the surface albedo available in the typical Lambertian surface albedo climatologies is more representative of the minimum value of the surface reflectivities that were observed (see, e.g., Lorente et al., 2018; Liu et al., 2020) – and it therefore underestimates the surface albedo for many of the scattering geometries. The disadvantage of the GLER approach is that it, at least for land-covered scenes, depends fully on the MODIS surface BRDF database. This limits the spectral usage to the seven wavelength bands of the MODIS BRDF product for land-covered scenes. For the retrieval of NO₂ and of cloud properties from the O₂–O₂ band, both performed in the spectral regime close to 466 nm, this is not a problem – but for many other retrievals it is.

A second example of a geometry-dependent surface LER database is the geometry-dependent effective Lambertian-equivalent reflectivity (GE_LER) database introduced in a recent paper by Loyola et al. (2020). The GE_LER approach does not depend on external data such as MODIS BRDF data, and it uses machine learning techniques to retrieve the surface reflectivity from level-1 data of the sensor (GOME-2, TROPOMI, or another UVN sensor). Like the GLER, the GE_LER provides daily maps of the surface properties. The GE_LER provides information for all surface types (land, ocean, snow/ice) in one database and covers the ultraviolet, visible, and near-infrared (UV–VIS–NIR) spectral region.

In this paper we introduce the directionally dependent Lambertian-equivalent reflectivity (DLER) of the Earth's surface derived from GOME-2 observations. The surface DLER is retrieved as a function of the viewing geometry and therefore describes the anisotropy of the surface reflectivity. The DLER approach is very different than the GLER approach in that we perform a retrieval directly on GOME-2 level-1 data, not relying on BRDF input (or any other input) from an external database. In this way the wavelength bands, 26 in total, can be chosen freely, allowing the resulting DLER database to support the retrieval of most atmospheric species. A difference compared to the GLER and GE_LER databases is that the directional dependence of the DLER is provided as a parameterisation of the viewing angle. It is not mapped on a satellite footprint and serves as a climatological dependence. The directional approach of the GOME-2 surface DLER is therefore applicable to all polar satellites with Equator crossing times close to that of GOME-2 (09:30 LT). This includes satellite instruments like GOME and SCIAMACHY, GOME-2 itself, and the future Sentinel-5 Ultraviolet, Visible, Near-infrared, and Shortwave infrared (UVNS) instrument scheduled for launch in 2023.

Like the GLER and GE_LER, the DLER is a Lambertian property and therefore can be used in situations where radiative transfer calculations include Lambertian surface reflection. The GOME-2 surface DLER database is an important improvement on the non-directional GOME-2 surface LER database that was described earlier (Tilstra et al., 2017). The transition from LER to DLER is the main topic of this paper including a study on the theoretical difference between DLER and BRDF data. Other improvements to the database are also described in this paper.

The paper is structured in the following way. Section 2 introduces the theory behind Lambertian-equivalent reflectivity (LER) and the new concept of directionally dependent LER (DLER). In Sect. 3 the theoretical difference between surface DLER and surface BRDF data is studied. Section 4 provides a short description of the GOME-2 instrument. The algorithm set-up, atmospheric correction, and the theoretical background of the improved surface DLER retrieval algorithm are described extensively in Sect. 5. Section 6 presents results and provides examples of the anisotropy of the Earth's surface according to the new GOME-2 surface DLER database. In Sect. 7 the DLER database is compared to MODIS BRDF data. Case studies and global comparisons are both performed, and the validation results are discussed. The paper ends with a summary and conclusions.

2 DLER theory

This section introduces the concept of a directionally dependent Lambertian-equivalent reflectivity (DLER) to describe the reflectivity of the Earth's surface. The following definition of the Earth reflectance is adopted in this paper:

$$\begin{array}{}\text{(1)}& R={\displaystyle \frac{\mathit{\pi }I}{{\mathit{\mu }}_{\mathrm{0}}E}}.\end{array}$$

In Eq. (1), the symbol I refers to the Earth radiance at the top of atmosphere (TOA; in ${\mathrm{Wm}}^{-\mathrm{2}}{\mathrm{sr}}^{-\mathrm{1}}{\mathrm{nm}}^{-\mathrm{1}}$). The symbol E refers to the incoming solar irradiance, perpendicular to the solar beam, at the TOA (given in ${\mathrm{Wm}}^{-\mathrm{2}}{\mathrm{nm}}^{-\mathrm{1}}$). The parameter μ₀ is a shorthand for μ₀=cos θ₀, with θ₀ the solar zenith angle. The shorthand for the viewing direction is μ=cos θ, with θ being the viewing zenith angle. The symbols for the viewing and solar azimuth angles are ϕ and ϕ₀, respectively.

2.1 Lambertian-equivalent reflectivity

The focus of this paper is on Lambertian surface reflection in combination with clear-sky atmospheric conditions. For these conditions, there exists a simple relationship between the Earth reflectance R and the (Lambertian) surface albedo A_s (Chandrasekhar, 1960):

$$\begin{array}{}\text{(2)}& R(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}},\mathit{\varphi },{\mathit{\varphi }}_{\mathrm{0}},A_{\mathrm{s}})=R^{\mathrm{0}}(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}},\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{0}})+{\displaystyle \frac{A_{\mathrm{s}}T(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}})}{\mathrm{1}-A_{\mathrm{s}}{s}^{\star }}}.\end{array}$$

In Eq. (2), the first term on the right is the so-called path reflectance R⁰. This is the atmospheric contribution to the Earth reflectance for a Rayleigh atmosphere which is bounded below by a non-reflecting surface. The second term in Eq. (2) is the surface contribution to the Earth reflectance. This term depends on the surface albedo A_s, on the total transmission T of the atmosphere, and on the spherical albedo s^⋆. The property s^⋆ is the spherical albedo of the Rayleigh atmosphere for illumination from below. The parameters R⁰, T, and s^⋆ can in principle be calculated using any radiative transfer model (see, e.g., Tilstra et al., 2012).

From a given measured reflectance R^obs, the surface albedo A_s can now be determined from Eq. (2):

$$\begin{array}{}\text{(3)}& A_{\mathrm{s}}={\displaystyle \frac{R_{\mathit{\lambda }}^{\text{obs}}-R_{\mathit{\lambda }}^{\mathrm{0}}}{T_{\mathit{\lambda }}(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}})+s_{\mathit{\lambda }}^{\star }(R_{\mathit{\lambda }}^{\text{obs}}-R_{\mathit{\lambda }}^{\mathrm{0}})}}.\end{array}$$

Both parameters $R_{\mathit{\lambda }}^{\text{obs}}$ and $R_{\mathit{\lambda }}^{\mathrm{0}}$ depend on μ, μ₀, and ϕ−ϕ₀ and so, in general, does A_s. When clear-sky conditions apply, the parameter A_s is the Lambertian-equivalent reflectivity (LER) of the surface.

Figure 1 (a) Illustration of the principle of Lambertian (isotropic) surface reflection. (b) Surface reflection distribution with a retroreflection lobe, representative of land surfaces covered by vegetation. In the DLER retrieval code, the orbit swath is divided into five viewing angle ranges, and for each segment the surface LER is determined in the usual way.

2.2 Directionally dependent surface LER

Traditional, non-directional surface LER databases are built on the assumption that all surface types act as Lambertian reflectors. That is, one assumes that the amount of light being reflected by the surface does not depend on the direction of incoming and reflected light. This is illustrated in Fig. 1a. The Lambertian assumption is, unfortunately, in many cases not justified. A more realistic description of the reflective properties of a surface requires a bi-directional reflectance distribution function (BRDF) (Nicodemus et al., 1992; Schaepman-Strub et al., 2006). The BRDF is a function of the incoming and outgoing directions. A hypothetical surface BRDF is shown in Fig. 1b. Here, the surface BRDF contains a retroreflection lobe resulting from increased reflection by vegetation in the backscattering direction.

In the retrieval algorithm of the traditional GOME-2 surface LER database (Tilstra et al., 2017), grid cells acts as storage containers in which all observations with fitting geolocation are stored, irrespective of viewing geometry and scene conditions. Statistical methods are then employed to identify the cloud-free scenes. For the directional GOME-2 surface DLER, the grid cell container is split into five sub-containers, each representing a certain viewing angle range (see Fig. 1b). The traditional retrieval algorithm is then run five times, deriving surface LER for each of the five viewing angle containers. The viewing angle dependence can then be analysed. This procedure is explained in Fig. 2.

Figure 2 Surface LER retrieved for a grid cell over the Sahara desert for the five viewing angle containers (indicated by the circles). The associated viewing angle θ_v is plotted on the horizontal axis. The vertical dotted lines indicate the centres of the viewing angle ranges. Colours indicate the selected wavelength bands. The curves are parabolic fits through the data points.

The coloured circles in Fig. 2 represent the surface LER retrieved by GOME-2 for a grid cell over the Sahara desert, for the five viewing angle containers and for the 26 wavelength bands defined in Sect. 5.1. The viewing angle θ_v presented on the horizontal axis is defined as follows.

$$\begin{array}{}\text{(4)}& {\mathit{\theta }}_{\mathrm{v}}=\left\{\begin{array}{ll}-\mathit{\theta },& \text{for the east viewing direction}\\ \mathit{\theta },& \text{for the west viewing direction}\end{array}\right.\end{array}$$

The centres of the viewing angle containers are indicated by the dotted vertical lines. Parabolic curves are fitted to the five retrieved surface LER values for all wavelength bands. Labels are provided for most of the wavelength bands.

From Fig. 2 it can be concluded that there is a clear dependence on θ_v. The parabolic fits suggest that the dependence may be parameterised as a function of the viewing angle θ_v in the following way:

$$\begin{array}{}\text{(5)}& A_{\text{DLER}}=A_{\text{LER}}+c_{\mathrm{0}}+c_{\mathrm{1}}\cdot {\mathit{\theta }}_{\mathrm{v}}+c_{\mathrm{2}}\cdot {\mathit{\theta }}_{\mathrm{v}}^{\mathrm{2}},\end{array}$$

where θ_v is negative on the east side of the orbit swath and positive on the west side of the orbit swath; see Eq. (4). The coefficients c₀, c₁, and c₂ are wavelength dependent and are calculated for each grid cell, provided that all five viewing angle segments are sufficiently filled with observations. For water bodies the coefficients c₀, c₁, and c₂ are set to zero because surface DLER and BRDF data cannot be cast into a climatology easily for water surfaces. This is because of the strong dependence on the viewing and solar angles for sun glint conditions and because of the dependence on parameters such as wind speed and chlorophyll concentration. The provided surface reflectance over water is therefore the standard minimum LER and more representative of the diffuse component of the surface reflection (Liu et al., 2020).

3 Theoretical study: DLER versus BRDF

In this section we study the theoretical difference between surface DLER and surface BRDF data. As explained in Sect. 1, BRDF and DLER are fundamentally different properties and as such cannot be expected to yield the same values or to take over the role of the other in radiative transfer calculations. Nevertheless, as the results in this section will show, for certain wavelength regimes BRDF and DLER data are numerically comparable. This allows for practical applications and validation of DLER by comparison with BRDF (and vice versa).

3.1 MODIS BRDF model

The MODIS Ross–Li surface BRDF model is a linear kernel-based BRDF model used to describe the surface reflectance of land surfaces. The surface anisotropy is described by two geometry-dependent kernels which have to be combined with the provided kernel coefficients if one wants to calculate the BRDF. The Li–Sparse kernel K_geo is the geometric kernel which describes the contribution of sunlit and shaded parts of a scene due to the presence of three-dimensional objects, typically trees. The Ross–Thick kernel K_vol is the volumetric kernel which describes the smaller-scale variation of the leaf canopy, i.e. the orientation of the leaves themselves.

The geometric and volumetric kernels are independent of wavelength. The wavelength dependence of the BRDF is contained entirely in the kernel coefficients. The expression for the surface reflectivity is as follows (Strahler et al., 1999).

$$\begin{array}{ll}{\displaystyle }{A}_{\mathrm{g}}(\mathit{\lambda },\mathit{\theta },{\mathit{\theta }}_{\mathrm{0}},\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{0}})=& {\displaystyle }{f}_{\text{iso}}\left(\mathit{\lambda }\right)\\ {\displaystyle }& {\displaystyle }+f_{\text{vol}}\left(\mathit{\lambda }\right)\cdot K_{\text{vol}}(\mathit{\theta },{\mathit{\theta }}_{\mathrm{0}},\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{0}})\\ \text{(6)}& {\displaystyle }& {\displaystyle }+f_{\text{geo}}\left(\mathit{\lambda }\right)\cdot K_{\text{geo}}(\mathit{\theta },{\mathit{\theta }}_{\mathrm{0}},\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{0}})\end{array}$$

The exact expressions needed to calculate K_vol and K_geo are provided in Appendix A. The coefficients f_iso, f_vol, and f_geo are the kernel coefficients of the isotropic, volumetric, and geometric contributions.

3.2 DLER and BRDF model calculations

For our model calculations we make use of the “Doubling-Adding KNMI” (DAK) radiative transfer code which will be described extensively in Sect. 5.3. Here we make use of surface reflection defined by a BRDF instead of Lambertian surface reflection as described by Lorente et al. (2018). We thereto provide DAK the three kernel coefficients (f_iso, f_vol, and f_geo) as defined in the MODIS ATBD (Strahler et al., 1999). Using this set-up, the TOA reflectance is calculated at a number of wavelengths, for the VZA and SZA nodes θ and θ₀ that are also part of the look-up tables (LUTs) described in Sect. 5.3, and for 360 equidistant values of the relative azimuth angle ϕ−ϕ₀. The surface elevation is set to zero (sea level) and the ozone column to 350 DU. As before, cloud and aerosols are not included. The calculations are performed monochromatically.

Next, the simulated surface DLER is retrieved from the simulated TOA reflectances using a similar set-up as the one described in Sect. 2. The only difference here is that the input reflectances are not measured but simulated by DAK, i.e., they are based on surface reflection described by the BRDF kernel coefficients (f_iso, f_vol, and f_geo). The differences between BRDF and DLER for all angles θ, θ₀, and ϕ−ϕ₀ are then analysed as a function of wavelength.

3.3 Analysis and discussion

The results for 772 nm are presented in Fig. 3 in the form of polar plots. The solar zenith angle was set to 32^∘. Figure 3a presents the BRDF, characterised by the kernel coefficients $(f_{\text{iso}},f_{\text{vol}},f_{\text{geo}})=(\mathrm{0.36},\mathrm{0.24},\mathrm{0.03})$. These kernel coefficients are representative of vegetated surfaces such as forests (Lorente et al., 2018). Figure 3b shows the TOA reflectance calculated by the DAK radiative transfer model (RTM). Note that the reflectance is similar to the BRDF. Figure 3c presents the retrieved DLER. The differences between BRDF and DLER appear to be small. This is confirmed by Fig. 3d, which presents the BRDF (red curve) and DLER (dotted blue curve) inside the principal plane ($\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{0}}=\mathrm{0}{}^{\circ }$ or 180^∘). Differences are found, but they are small even for large viewing zenith angles.

Figure 3 (a, b) Surface BRDF at 772 nm for a solar zenith angle of 32^∘ and the resulting simulated TOA reflectance. The BRDF kernel coefficients (f_iso, f_vol, and f_geo) were set to (0.36, 0.24, and 0.03), representative of vegetation. (c, d) Retrieved surface DLER and BRDF–DLER in the principal plane. In the principal plane, where ϕ−ϕ₀ is 0^∘ or 180^∘, exact backscattering occurs at a viewing angle of 32^∘.

These results are not unexpected because at 772 nm the Rayleigh optical thickness is quite low (about 0.02), so scattering in the atmosphere is relatively weak. This means that (single) surface reflection dominates and that the reflectance at the TOA is similar to the BRDF of the surface. Moreover, in this situation the retrieved surface DLER and BRDF are similar. Note that the behaviour of the BRDF for extreme viewing zenith angles in the forward scattering direction is suspicious because the BRDF becomes negative for viewing zenith angles close to 90^∘. In the DAK RTM, the surface BRDF is therefore not allowed to become negative.

Figure 4 Similar to Fig. 3 but now for 463 nm. Note that the BRDF kernel coefficients are different from the ones in Fig. 3. The relative differences between DLER and BRDF are now larger, especially near the hot spot and for large viewing zenith angles.

The situation changes quite a bit at 463 nm (see Fig. 4). Figure 4 shows that at this wavelength the typical value of the BRDF is much lower than at 772 nm. Because of increased Rayleigh scattering in the atmosphere, the TOA reflectance is now very different from the surface BRDF. The retrieved DLER shows quite some differences compared to the BRDF. This is caused by an increased occurrence of scattering in the atmosphere (Rayleigh optical thickness ∼0.2 at 463 nm) and multiple scattering via the surface. The effects may seem to be modest in an absolute sense, but relative to the typical value of the BRDF (∼0.03) they can be quite large, depending on the viewing and solar angles that are involved. For instance, inside the “hot spot” the differences are 0.003, of which the BRDF is 0.045 and the DLER is 0.042. The differences therefore can go up to 6 % at this wavelength.

Based on the results we can distinguish three wavelengths regimes. For λ>1000 nm the functional behaviour of DLER and BRDF is nearly identical and interchangeable. For $\mathrm{500}<\mathit{\lambda }<\mathrm{1000}$ nm the DLER and BRDF values are similar, and DLER and BRDF can be interchanged in most practical situations. For example, at 555 nm the difference in the hot spot region is less than 3 % for ${\mathit{\theta }}_{\mathrm{0}}=\mathrm{32}{}^{\circ }$ and less than 8 % for ${\mathit{\theta }}_{\mathrm{0}}=\mathrm{60}{}^{\circ }$. For λ<500 nm, however, DLER and BRDF differ by too much for them to take over each other's role. For example, in the UV at 380 nm the differences go up to 12 % for ${\mathit{\theta }}_{\mathrm{0}}=\mathrm{32}{}^{\circ }$ and even 30 % for ${\mathit{\theta }}_{\mathrm{0}}=\mathrm{60}{}^{\circ }$. It should be noted that the absolute differences in these cases are small (maximum of the order of 0.01). Please note that vegetated surfaces have a relatively strong surface anisotropy compared to most other surface types such as desert. The provided numbers therefore represent worst-case situations.

Depending on the application, using BRDF instead of DLER, or DLER instead of BRDF, can be acceptable even below 500 nm. For the validation study presented in Sect. 7 we will, however, restrict ourselves to λ>500 nm.

4 Description of GOME-2

GOME-2 (Munro et al., 2016) is the successor of the Global Ozone Monitoring Experiment (GOME) (Burrows et al., 1999). It is a remote sensing spectrometer that measures the Earth's radiance and the solar irradiance, covering the wavelength range between 240 and 800 nm. The spectral resolution varies between 0.3 and 0.5 nm. Like its predecessor GOME, GOME-2 performs scans of the Earth in a motion from east to west in 4.5 s (forward scan) and back from west to east in 1.5 s (backscan). This motion is achieved via the rotation of an internal scanner mirror. The orbit swath is 1920 km wide and the measurement footprint in the forward scan is 80 km×40 km (across track×along track). In about 1.5 d every location on the Earth's surface is observed.

Next to the spectral measurements of radiance and solar irradiance, also the polarisation of the light is measured. This is done by on-board polarisation measurement devices (PMDs) which measure the state of atmospheric polarisation in 15 wavelength bands. The polarisation information is subsequently used to perform a correction for polarisation on the detected signals. Note that GOME-2 is sensitive to the polarisation of the incoming light since it is not equipped with a polarisation scrambler. The information from the PMD bands was used by us to also derive a surface DLER database based on the PMD bands. This PMD-based database is not described explicitly in this paper.

The first GOME-2 instrument was launched on 19 October 2006 as part of the MetOp-A satellite platform. Identical versions of the first GOME-2 instrument were launched on board the MetOp-B and MetOp-C satellites, with launch dates of 17 September 2012 and 7 November 2018, respectively. All three MetOp satellites were put into near-polar, Sun-synchronous orbits at an altitude of 820 km and with an orbital period of about 101 min. The local Equator crossing time is 09:30 LT for the descending node for all three satellite platforms but with different phasing. The MetOp series of satellites is expected to continue operations beyond the year 2027.¹

The GOME-2 instruments were designed to perform global observations of trace gases for environmental and meteorological applications and climate monitoring. Trace gases that are retrieved are ozone, NO₂, BrO, SO₂, HCHO, OClO, and H₂O. Next to trace gases also cloud, aerosol, and surface properties are retrieved. A complete overview of the available GOME-2 products is presented in Hassinen et al. (2016).

5 Algorithm set-up and atmospheric correction

The DLER algorithm set-up is in many aspects similar to the LER algorithm set-up described in Tilstra et al. (2017). That is, the Earth reflectance spectrum is transformed into a number of reflectance bands, which are converted into scene LER values by applying the atmospheric correction outlined in Sect. 2.1. After these steps, the observed scene LER values of a specific month (but from all available years) are distributed onto a latitude and longitude grid which represents the Earth's surface for that specific calendar month. In this step, observations containing absorbing aerosols are filtered out using the Absorbing Aerosol Index (AAI). For each grid cell the distribution of scene LER values is then analysed statistically to find the cloud-free observations. This is done in two ways. In the first method the so-called MIN-LER is retrieved, which is the 1 % cumulative value of the scene LER distribution. The second method retrieves the so-called MODE-LER field. The MODE-LER is found from the mode of the scene LER distribution, which is a well-defined maximum for arid (desert) surfaces and snow/ice surfaces. For the other surface types, the mode cannot be used, and the MIN-LER result is copied. Both MIN-LER and MODE-LER fields are present in the surface LER database. After these steps, post-processing corrections are performed that remove issues such as gaps and residual cloud contamination. The above steps and procedures have been described extensively in Tilstra et al. (2017). There are, however, a number of important improvements and extensions in the current algorithm:

1. The list of wavelength bands was extended with wavelength bands at 328, 585, 685, 697, and 712 nm. The wavelength bands at 685, 697, and 712 nm were introduced specifically to support cloud and aerosol retrieval near the O₂-B band, as explained in Desmons et al. (2019). See Sect. 5.1 and Table 1 for details.

2. Spectral calculations were introduced for some of the wavelength bands, and absorption by oxygen and water vapour was included in the way described in Sect. 5.2.

3. The spatial resolution of the database fields was increased using dynamic gridding. This dynamic gridding approach is explained in Appendix B.

4. Data from both MetOp-A and MetOp-B were used from the period 2007–2018, covering more than 10 years of observations. The larger amount of data used is beneficial for the quality of the climatology. Data from MetOp-A were used only until 2013 because in July 2013 the GOME-2A orbit swath was reduced from the standard 1920 to 960 km. The reduction of the viewing angle range would have impacted the non-directional LER since it would then have been biased towards the LER values of the inner part of the orbit swath.

5. The database now offers directionally dependent surface LER (DLER). This means that the anisotropy of the surface reflection, often called the BRDF effect, is contained in (and described by) the DLER database. The provided DLER is an approximation in the sense that the second-order polynomial approach presented in Sect. 2.2 in combination with the five angular bins of about 20^∘ each will not be able to catch the angular variability in the DLER for all surface types and situations. In particular, for vegetated surfaces the hot spot will not in all circumstances and geometries be represented well. Also, the DLER database is in principle representative only of the geometry of the GOME-2 orbit.

5.1 Selection of wavelength bands

Table 1 provides a list of the chosen wavelength bands, as well as their central wavelength and bandwidth. Note that the wavelength bands at 328, 585, 685, 697, and 712 nm were not present in the previous version of the GOME-2 surface LER database (Tilstra et al., 2017, their Table 2). Most of the wavelength bands are 1 nm wide. The wavelength bands are therefore narrow enough to be considered monochromatic but also wide enough to effectively minimise the impact of the Ring effect (Chance and Spurr, 1997). The reflectances for the wavelength bands are calculated from the reflectances measured by the individual detector pixels that fall within the wavelength band. A boxcar weighting function w is applied to each detector pixel reflectance. This weighting function is defined as follows.

$$\begin{array}{}\text{(7)}& w_i^j=\left\{\begin{array}{ll}\mathrm{1},& \text{for}|{\mathit{\lambda }}_i-{\mathit{\lambda }}_j^{\mathrm{c}}|\le {\mathit{\omega }}_j\\ \mathrm{0},& \text{for}|{\mathit{\lambda }}_i-{\mathit{\lambda }}_j^{\mathrm{c}}|>{\mathit{\omega }}_j\end{array}\right.\end{array}$$

In this equation, λ_i is the wavelength of detector pixel i, ${\mathit{\lambda }}_j^{\mathrm{c}}$ is the central wavelength of wavelength band j, and 2ω_j is the width of wavelength band j. Normalisation of the resulting band reflectance is performed by dividing the result with the number of participating detector pixels, denoted by N_j.

Table 1 Definition of the wavelength bands and details of the radiative transfer calculations for atmospheric correction.

The reflectance calculations are performed using spectral band integration or monochromatically. For all wavelength bands absorption by ozone, NO₂, and O₂–O₂ is included. Absorption by oxygen and/or water vapour is included for only some of the wavelength bands.

Most of the wavelength bands are positioned in the continuum parts of the spectrum, avoiding absorption bands as much as possible. This is essential because having to take absorption by atmospheric species into account complicates the radiative transfer calculations considerably. For wavelength bands located in the continuum monochromatic simulations are sufficient. This is not the case for a number of wavelength bands which are affected too much from absorption by trace gases. These wavelength bands (see Table 1) require a spectral handling of the radiative transfer calculations. This approach is described in Sect. 5.2.

5.2 Absorption by trace gases

Three examples of the impact of absorption by oxygen, water vapour, and ozone are given in Fig. 5. Figure 5a shows the situation for the wavelength band at 758 nm, positioned just in front of the O₂-A band, while Fig. 5b shows the situation for the wavelength band near 697 nm, which is spectrally surrounded by water vapour absorption lines. Figure 5c presents the situation for the wavelength band at 328 nm, where ozone absorption is quite variable over the extent of the wavelength band. The black curves represent the simulated reflectance spectra. These spectra were calculated for clear-sky conditions, for a surface albedo A_s=0.5 at sea level, for nadir view and local noon ($\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}=\mathrm{0}{}^{\circ }$), for an ozone column of 350 DU, and for a water vapour column of 4.0 g cm⁻². For comparison, the horizontal green curves represent the reflectance spectra without taking absorption by oxygen and/or water vapour into account in the radiative transfer calculations.

Figure 5 Simulated reflectance spectra (in black) relevant to the wavelength bands at 758, 697, and 328 nm. Panel (a) shows the very small impact of oxygen absorption near 758 nm, while (b) shows the larger impact of water vapour absorption around 697 nm. Panel (c) presents the situation of ozone absorption near 328 nm. For comparison, the horizontal green curves indicate reflectance spectra simulated without absorption by oxygen and water vapour. The vertical green lines indicate the spectral positions of the detector pixels that make up the wavelength bands (see Table 1). The blue curves indicate the response functions of the wavelength bands based on the indicated detector pixels and their individual slit functions.

For the wavelength band at 758 nm the impact of oxygen absorption is obviously very small. On the other hand, for the wavelength band near 697 nm the impact of water vapour absorption is much larger. A monochromatic calculation is clearly not sufficient in this case. To proceed, we first define G_j(λ), the spectral response function of wavelength band j, as a weighted superposition of the slit functions of the individual detector pixels by using the boxcar weighting defined in Eq. (7). That is, for the response function G_j we have

$$\begin{array}{}\text{(8)}& G_j\left(\mathit{\lambda }\right)={\displaystyle \frac{\mathrm{1}}{N_j}}\sum _{i=\mathrm{1}}{w}_i^j\cdot S_i\left(\mathit{\lambda }\right),\end{array}$$

where S_i is the normalised slit function of detector pixel i from the appropriate spectral band and N_j the number of detector pixels that make up wavelength band j. The resulting response functions G₇₅₈, G₆₉₇, and G₃₂₈ are presented in Fig. 5 as blue curves, in arbitrary units. The vertical green lines indicate the wavelengths λ_i of the detector pixels that contribute to the reflectance of the wavelength band. Next, we calculate the simulated band reflectance. For this we first need to simulate the spectrum surrounding the wavelength band at a high spectral resolution. We use a spectral resolution of 0.01 nm. The spectral sampling is then increased by a factor of 100 using Akima interpolation (Akima, 1970). This allows for accurate numerical integration, and the resulting expression for the simulated band reflectance is

$$\begin{array}{}\text{(9)}& R_j^{\text{sim}}=\sum _k\mathrm{\Delta }\mathit{\lambda }\cdot G_j\left({\mathit{\lambda }}_k\right)\cdot R^{\text{sim}}\left({\mathit{\lambda }}_k\right),\end{array}$$

where the summation over k involves a summation over the wavelengths λ_k and $\mathrm{\Delta }\mathit{\lambda }={\mathrm{10}}^{-\mathrm{4}}$ nm.

The impact of neglecting absorption by oxygen and/or water vapour and using monochromatic calculations can now be calculated. For the 758 nm case given in Fig. 5 this effect is only 0.003 on the reflectance and about the same for the surface LER. This wavelength band could therefore be treated monochromatically. For the 697 nm case, however, the effect is 0.018, which is too high to justify monochromatic calculations. For the 328 nm wavelength band, adopting monochromatic calculation would lead to an error of −0.027. For the wavelength bands at 328, 697, 712, 758, and 772 nm we use Eq. (9) to calculate the simulated band reflectance. In Table 1 this is indicated by the label “S” in the fifth row. For the other wavelength bands we use monochromatic calculations, indicated by “M” in the fifth row of Table 1. For the 758 and 772 nm wavelength bands a monochromatic calculation would have sufficed, but because of their strong importance to cloud and aerosol retrieval using the O₂-A band, we decided to go further than necessary by adopting spectral calculations.

5.3 Radiative transfer calculations and LUTs

For the radiative transfer calculations we make use of the radiative transfer code “Doubling-Adding KNMI” (DAK) (de Haan, 1987; Stammes, 2001). The DAK code is able to calculate all four components of the Stokes vector (van de Hulst, 1981; Hovenier et al., 2004), and in its minimal set-up it features molecular scattering and Lambertian surface reflection, but the user can decide to include many other features such as scattering by clouds and/or aerosols, absorption by various trace gases, and surface reflection defined by a BRDF. The extension to BRDF was described by Lorente et al. (2018). In the calculations we did not include clouds and aerosols and adopted Lambertian surface reflection. Polarisation is included in the calculations. We used a standard mid-latitude summer (MLS) atmosphere (Anderson et al., 1986) for the atmospheric profiles and included absorption by ozone, NO₂, and O₂–O₂ for all wavelength bands. For some of the wavelength bands we also included absorption by oxygen and water vapour (see Table 1).

For all 26 wavelength bands look-up tables (LUTs) were created. The LUTs were made for 7 ozone column values (50, 200, 300, 350, 400, 500, and 650 DU), for 10 surface heights (ranging from 0 to 9 km in steps of 1 km), for water vapour columns of 0 and 4 g cm⁻², and for 42 non-equidistant values of μ and μ₀. The dependence on the relative azimuth angle ϕ−ϕ₀ can be treated analytically. To explain, because the simulations represent clear-sky Rayleigh atmospheres, the Fourier expansion of the reflectance in terms of the relative azimuth angle ϕ−ϕ₀ ends after only three terms. More specifically, we have

$$\begin{array}{}\text{(10)}& R^{\mathrm{0}}=a_{\mathrm{0}}(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}})+\sum _{i=\mathrm{1}}^{\mathrm{2}}\mathrm{2}{a}_i(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}})\cos i(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{0}}).\end{array}$$

We therefore do not store the reflectances R⁰ in the LUTs but instead store the Fourier coefficients a₀, a₁, and a₂. The reflectance R⁰ can be calculated from these. The LUTs contain the parameters a₀, a₁, a₂, T, and s^⋆.

6 Results

6.1 Surface anisotropy

Examples of the magnitude of the DLER surface anisotropy in GOME-2 data are provided in Fig. 6. In the left column the traditional non-directional GOME-2 surface LER is shown. The parameter presented in the right column is the surface anisotropy parameter, defined here as the difference between the GOME-2 surface DLER at viewing angles θ_v of +45^∘ (west viewing direction) and −45^∘ (east viewing direction). For the GOME-2 orbit this parameter is a good indicator of the magnitude and range of the surface anisotropy in the GOME-2 orbit swath. The results in Fig. 6 are presented for calendar month March and for the wavelength bands at 772, 670, and 555 nm. At 772 nm the surface anisotropy parameter can be as large as 0.2 for vegetated areas. For the typical desert areas the differences are much smaller (0.05–0.10) because non-vegetated surfaces are usually more isotropic than vegetated areas. The surface anisotropy parameter for snow/ice surfaces has the opposite sign. The values over the vegetated areas correspond to percentages of 50 %–125 %, in agreement with what was found already by Lorente et al. (2018). The magnitude of the surface anisotropy which is present in the GOME-2 surface DLER is therefore quite substantial.

Figure 6 (a, c, e) Global maps of the GOME-2 surface LER for calendar month March and for 772, 670, and 555 nm. (b, d, f) Global maps of the surface anisotropy parameter, defined as the difference between GOME-2 surface DLER at viewing angles of +45^∘ and −45^∘. The surface anisotropy can be large, especially for vegetated surfaces at wavelengths beyond 700 nm. Over the oceans only non-directional surface LER is provided, as explained in Sect. 2.2.

At 670 and 555 nm the surface anisotropy parameter over vegetation is much lower than at 772 nm. However, this is mainly caused by the fact that the surface reflectance at these wavelength bands is also much lower than at 772 nm. The percentages are more or less the same, in the range of 50 %–125 %. For desert areas the anisotropy parameter is slightly smaller at 670 and 555 nm than at 772 nm. However, the percentages are similar for all three wavelength bands, about 10 %–20 %. For snow/ice surfaces the anisotropy parameter is not depending much on the wavelength. For 555 nm the values are slightly smaller. This is probably because of increased Rayleigh scattering, which results in a more diffuse illumination of the surface. The surface anisotropy parameter varies mostly between −0.1 and −0.3 for snow/ice surfaces.

6.2 Dependence on surface type and time

The directional dependence of the surface DLER was studied for the nine surface types defined in Table 2. All nine surface types represent land surfaces. Constraints were set on latitude and longitude and, more importantly, on surface type using the Matthews land usage database (Matthews, 1983). Furthermore, coastal areas were excluded and so were all grid cells that contained snow/ice, except for the “Antarctica” and “Greenland” surface types. The results are summarised in Fig. 7, which presents the GOME-2 surface DLER as a function of the viewing angle θ_v defined in Eq. (4). The solid coloured curves represent the surface DLER at 772 nm, averaged over the grid cells of the surface type region as defined in Table 2. The legend in the “Antarctica” window explains to which months the curves belong. The grey curves in Fig. 7 are there to provide an indication of the spread in the surface DLER. The spread is defined as 2.35 times the standard deviation in the data. The legend in the “Greenland” window explains to which months the grey curves belong.

Table 2 Definition of the surface type regions studied in Fig. 7. The symbol “–” indicates that no constraint was set on the longitude.

Figure 7 Surface DLER at 772 nm versus viewing angle for nine surface types and four calendar months. The coloured curves represent the average GOME-2 surface DLER. The grey curves provide an indication of the spread in surface DLER over the selected regions.

The three desert surface types (“Sahara desert”, “Arabian Peninsula”, and “Australian desert”) show a very similar dependence on the viewing angle. The overall dependence agrees well with that of the single grid cell shown in Fig. 2. The Australian desert deviates slightly from the other two desert regions because of the much lower values and the larger variability with respect to the calendar month. The latter observation may be partly explained by the different solar zenith angles, but more likely it is caused by the fact that the Australian desert contains more vegetation than the other two desert surface types. The snow/ice surface types (“Antarctica” and “Greenland”) show that the highest value for the surface DLER is reached for the eastward looking direction and not for the westward looking direction as for the other surface types presented in Fig. 7. There is, at least for Antarctica, a mild dependence on the calendar month. Note that for certain months the surface DLER results are not plotted for these regions. This is because the regions are covered in polar night during these months. The surface DLER is available for these months, but it is a replacement based on other months, which is why the results are not plotted in Fig. 7.

The four remaining surface type regions (“Amazonian tropical rainforests”, “Asian (sub-)tropical forests”, “Deciduous forests”, and “Grasslands”) show a large dependence on viewing angle. For example, for the month of May (brown curve) the average surface DLER in the Amazonian region varies from 0.21 for θ_v of −55^∘ to 0.36 for θ_v of +55^∘. For the month of November (blue curve) the increase is from 0.22 to 0.49, which is more than a factor of 2. The variability in time is the largest at the westward looking viewing direction.

The Asian (sub-)tropical forests and the Deciduous forests on the other hand show a temporal variability which is similar for the entire viewing angle range. Grasslands show a low temporal variability. For all four vegetated surfaces the anisotropy of the surface reflection is large. For these surface types the advantage of using DLER instead of LER is therefore substantial, but also for desert and snow/ice surfaces there is a significant improvement. Results for other wavelength bands can be found in Figs. S1–S3 in the Supplement.

7 Validation and discussion

In this section the GOME-2 surface DLER database is compared to the established MODIS surface BRDF product. This is done in two ways. First, in Sect. 7.1, case studies will be performed to analyse the directional behaviour of the two surface reflectivity products. Then, in Sect. 7.2, global comparisons will be presented. In both sections we make use of the MODIS MCD43C2 snow-free product, which provides surface BRDF for snow-free land scenes for seven of the MODIS bands. We select MODIS band 1, centred around 645 nm, as a reference for the 640 nm wavelength band of the GOME-2 surface DLER database. The choice for MODIS band 1 is based on the fact that (i) it is close enough to one of the DLER wavelength bands, and (ii) based on the results presented in Sect. 3 we may expect only small differences between DLER and BRDF for wavelengths longer than 600 nm.

7.1 Case studies

In Fig. 8 we present the results from a comparison between GOME-2 surface DLER and MODIS surface BRDF for three surface type cases. The three reference sites (Amazonian rainforest; equatorial Africa; Libyan desert) were selected primarily on the basis of their homogeneity. Homogeneity is important because the MODIS MCD43C2 product is provided at a 5×5 times higher spatial resolution than the GOME-2 surface DLER database. In all three cases we selected a one-by-one degree latitude and longitude box, containing 16 grid cells from the GOME-2 surface LER database and 400 grid cells from the MODIS MCD43C2 database. The selected grid cells supply all the surface reflectivity parameters that are needed to calculate DLER and BRDF from the equations that were introduced earlier.

Figure 8 Surface reflectivity around 640 nm according to the GOME-2 surface DLER database (black curves) and the MODIS surface BRDF product (blue curves), as a function of the viewing angle, for the Amazonian rainforest (a), equatorial Africa (b), and the Libyan desert (c). Results are representative of 15 March 2008 and correspond to a one-by-one degree latitude and longitude box with the indicated central coordinates.

We then feed the geometry-dependent DLER and BRDF equations with artificial but realistic GOME-2 viewing and solar angles. The viewing angle θ_v is varied between −55^∘ and +55^∘ to simulate the scanning motion of the GOME-2 instrument. This is already sufficient information to calculate the DLER for the 16 selected DLER grid cells (see Eq. 5). For the BRDF, the viewing zenith angle is then automatically known ($\mathit{\theta }=\left|{\mathit{\theta }}_{\mathrm{v}}\right|$; see Eq. 4), but the solar zenith angle θ₀ and relative azimuth angle ϕ−ϕ₀ also need to be known. The solar angles θ₀ and ϕ₀ can be determined from solar position calculations (Michalsky, 1988) by specifying an hour angle determined from the GOME-2 Equator overpass time of 09:30 LT, taking into account the change in hour angle because of the displacement in the latitude direction (i.e., caused by the rotation of the Earth) and in longitude direction (because of the scanning motion from east to west). The viewing azimuth angle ϕ is quite a constant factor (apart from a 180^∘ jump when GOME-2 scans past the exact nadir direction) and was determined from GOME-2 data. With all artificial angles known, the MODIS kernels can be calculated, and subsequently the surface BRDF can be calculated using Eq. (6) for all 400 MODIS MCD43C2 grid cells.

In Fig. 8a the scene that is studied is located over the Amazonian rainforest. The blue curves represent the surface BRDF from MODIS band 1 from the 400 MODIS BRDF grid cells. Note that even for what we consider homogeneous scenes there is already quite some variability between the grid cells. The GOME-2 surface DLER at 640 nm is presented in black. The agreement is good, both qualitatively in terms of the directional dependence and in absolute sense. Figure 8b presents the case of a scene in equatorial Africa. Here the agreement is again good with the correct dependence on the viewing angle. For the east viewing directions (${\mathit{\theta }}_{\mathrm{v}}<-\mathrm{50}{}^{\circ }$), however, the agreement seems to be a little less good. It should be noted that there is considerable variability in the 400 blue curves and that some of the blue curves also show the same upward bend for ${\mathit{\theta }}_{\mathrm{v}}<-\mathrm{50}{}^{\circ }$ as the black curves.

Finally, in Fig. 8c the scene studied is over the Libyan desert. The Libyan desert is known to be rather stable and is often used as a calibration reference site (e.g. Tilstra et al., 2005). The variability in the approximately 100 km by 100 km large box is low. The agreement between DLER and BRDF is good, with small 3 %–4 % differences for the east viewing directions (${\mathit{\theta }}_{\mathrm{v}}<-\mathrm{50}{}^{\circ }$). The lower performance at the most extreme viewing angles could be a result of the fact that in the parabolic fitting procedure explained in Sect. 2.2 inaccuracies are not corrected at the swath ends, but they are in the centre of the swath. Also, the parameterisation used for the DLER is a second-order polynomial, so higher-order dependencies are not described well. Another explanation could be background aerosol scattering. There is no explicit filtering or correction for aerosol scattering in the retrieval code. Aerosol scattering would have the largest impact at the extreme viewing angles.

We conclude that the DLER follows the correct directional behaviour. Results for other wavelength bands can be found in Figs. S4 and S5 in the Supplement. In the next section global comparison are performed to be able to draw more quantitative conclusions.

7.2 Global comparisons

This section presents results from global comparisons between GOME-2 surface DLER and MODIS surface BRDF. This time, the scattering geometry defined by θ, θ₀, and ϕ−ϕ₀ is prescribed by real GOME-2 observations. For these GOME-2 observations the surface DLER is calculated from the closest grid cell of the GOME-2 surface DLER database for the appropriate viewing angle θ_v (see Sect. 7.1). The MODIS surface BRDF is taken from the MODIS MCD43C2 snow-free product and based on the exact 25 grid cells that coincide with the grid cell of the GOME-2 surface DLER database. The MODIS surface BRDF is calculated for the prescribed θ, θ₀, and ϕ−ϕ₀ and averaged over the 25 grid cells. Only observations from the descending forward scan and over land are accepted. Scenes over ocean or scenes containing snow or ice are skipped. Observations with absolute latitudes above 60^∘ are also skipped.

Figure 9 Global comparisons between GOME-2 surface (D)LER and MODIS surface BRDF for the period 10–19 May 2019. (a, b) GOME-2 non-directional surface LER at 640 nm versus MODIS surface BRDF from band 1 (centred around 645 nm) for eastern and western sides of the orbit swath (see main text). (c, d) GOME-2 directional surface DLER versus MODIS surface BRDF. For the western side of the orbit swath the directional database agrees much better with MODIS than the non-directional one.

A typical outcome for 640 nm is presented in Fig. 9. The result was obtained for 10 d of GOME-2B observation geometries (10–19 May 2019). In Fig. 9a and b the traditional, non-directional GOME-2 surface LER is presented against the MODIS surface BRDF for the eastern side of the orbit swath (IndexInScan 1–8; ${\mathit{\theta }}_{\mathrm{v}}<-\mathrm{23}{}^{\circ }$; Fig. 9a) and the western side of the orbit swath (IndexInScan 17–24; ${\mathit{\theta }}_{\mathrm{v}}>+\mathrm{23}{}^{\circ }$; Fig. 9b). For the eastern side of the orbit swath the correlation between surface LER and MODIS BRDF is quite fair. The linear fit has a slope of 0.845, which admittedly deviates from 1, but Pearson's correlation coefficient r is 0.955, indicating good correlation. The standard deviation of the data points with respect to the linear fit, σ, amounts to 0.027. This value of σ is in line with uncertainties of ∼0.02 that were reported for the GOME-2 surface LER database and for a few other surface LER databases (Tilstra et al., 2017).

For the western side of the orbit swath, Pearson's r still suggests a reasonable correlation between the two datasets, but the linear fit deviates quite a bit more from the one-to-one relationship. Also, the scatter plot suggests that there is no pure linear relationship between the two databases. The explanation for this behaviour is the bias in the traditional LER databases towards the geometries with the lowest surface LER values. For GOME-2 this affects mostly the western side of the orbit swath, as reported earlier by Lorente et al. (2018).

In Fig. 9c and d the results are presented for the comparison between the directional GOME-2 surface DLER and MODIS surface BRDF. For the eastern side of the orbit swath (Fig. 9c) there are no clear changes compared to the situation in Fig. 9a. The slope and intercept have improved only marginally. For the western side of the orbit swath (Fig. 9d), however, the correlation has improved considerably. Both slope and intercept have improved, Pearson's r indicates higher correlation, and σ went down mildly. More importantly, the eastern and western side of the orbit swath now appear to perform equally well.

The analysis presented in Fig. 9 was repeated for each calendar month for the period 2012–2019 to search for time dependencies in the results. Clear time dependencies, either seasonal or annual, were not found. We also performed the analysis using the SCIAMACHY surface LER database instead of the GOME-2 surface DLER database. The SCIAMACHY surface LER database is by definition a non-directional database, so we could only produce results as in the top row of Fig. 9. The resulting scatter plots were very comparable. For instance, the linear fit to the data had a slope of 0.860 and an intercept of 0.005, which agrees with the numbers in Fig. 9 (0.845 and 0.005, respectively). This may suggest that the deviation of the slope from 1 is not just specific for GOME-2 but specific for the differences between the LER databases in general and MODIS BRDF.

In the past, small radiometric calibration errors in the GOME-2 level-1 data have been reported (Cai et al., 2012; Wu et al., 2014; Tilstra et al., 2014, 2017). However, these were mainly found for the UV wavelength range and not so much for the particular wavelength that was studied here (640 nm). Also, other studies have reported good agreement with Advanced Very High Resolution Radiometer (AVHHR) and Advanced Along-Track Scanning Radiometer (AATSR) for the wavelength range 630–670 nm (Latter et al., 2011). This suggests that the results and conclusions of this section were not significantly influenced by calibration errors in the GOME-2 data.

Results for other wavelength bands can be found in Figs. S6 and S7 in the Supplement.

8 Conclusions

In this paper we introduced the directionally dependent Lambertian-equivalent reflectivity (DLER) of the Earth's surface, retrieved from GOME-2 observations. This directional GOME-2 surface DLER database is a major update of the previous non-directional GOME-2 surface LER database in the sense that it describes the anisotropy of surface reflection, while the traditional LER database considers surface reflection to be isotropic. The DLER database can be used in atmospheric trace gas, aerosol, and cloud retrieval algorithms, just like the previous LER database. The retrieval of DLER was described and the anisotropy of the surface reflection could be studied. The anisotropy is especially large for the longer wavelengths and for vegetated surfaces.

Other improvements to the GOME-2 surface DLER database were also described. These include additional wavelength bands, an improved atmospheric correction taking absorption by oxygen and/or water vapour into account, a higher quality due to the use of more mission data, and a higher spatial resolution. The higher spatial resolution was achieved without compromising the quality by adopting a dynamic gridding approach. However, the main improvement is in the directional nature of the DLER.

To analyse the newly defined property, we conducted a series of radiative transfer simulations to study the theoretical differences between DLER and BRDF. The study showed that DLER and BRDF are different for the shorter wavelengths (λ<500 nm). Here the DLER is meant to be used only in combination with a radiative transfer code that includes Lambertian surface reflection. Lambertian surface reflection is the simplest form of surface reflection and probably the most used in practice. For the longer wavelengths, the differences between DLER and BRDF are small ($\mathrm{500}\mathrm{nm}<\mathit{\lambda }<\mathrm{1000}$ nm) to negligible (λ>1000 nm), and DLER can effectively be used as a BRDF (and BRDF as DLER).

This conclusion allowed us to compare the GOME-2 surface DLER at 640 nm with MODIS surface BRDF from MODIS band 1 (centred around 645 nm). A few case studies illustrate that the angular dependencies of DLER and BRDF are indeed comparable and that there is good agreement between DLER and BRDF. After that, extensive global comparisons confirm that there is indeed good systematic agreement. Moreover, the comparison with MODIS BRDF is also performed using the traditional, non-directional LER database. This LER database performs rather badly at the western side of the GOME-2 orbit swath. This is in line with findings by Lorente et al. (2018), who found that traditional, non-directional LER databases underestimate the surface reflection for certain viewing geometries. In particular, for the GOME-2 orbit geometry it was found that the surface reflection at 772 nm was underestimated by a factor of 2 at the western side of the orbit swath. This is in agreement with our findings. For the western side of the orbit swath, the DLER performs considerably better than the LER database.

The directional DLER database is, in summary, an important improvement on the traditional non-directional LER databases that are often used in atmospheric retrieval applications. The GOME-2 surface DLER database can be used as an input parameter for atmospheric retrieval algorithms working on data from all polar satellites with an Equator crossing time close to that of GOME-2. This includes instruments such as GOME, SCIAMACHY, and GOME-2 itself, as well as the future Sentinel-5 UVNS instrument which is scheduled for launch in 2023.

Appendix A: Kernels for the Ross–Li BRDF model

This appendix lists the equations needed to calculate the kernels that make up the Ross–Li BRDF model of surface reflectance. Proper derivations of the Ross–Thick and Li–Sparse kernels can be found in Wanner et al. (1995).

A1 Ross–Thick volumetric kernel

The Ross–Thick volumetric scattering kernel is defined in the following way (Roujean et al., 1992):

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

In Eq. (A1), θ refers to the viewing zenith angle and θ₀ to the solar zenith angle. The angle ξ is defined according to

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

where ψ and ψ^′ are the viewing and solar azimuth angles following the definition in Strahler et al. (1999). Exact backscattering ($\mathit{\xi }=\mathrm{0}{}^{\circ }$) occurs for $\mathit{\psi }-{\mathit{\psi }}^{\prime }=\mathrm{0}{}^{\circ }$, which agrees with the definition used for the GOME-2 data products.

A2 Li–Sparse geometric kernel

The Li–Sparse geometric scattering kernel (Li and Strahler, 1986) is defined as

$$\begin{array}{}\text{(A3)}& K_{\text{geo}}=O-\sec {\mathit{\theta }}^{\star }-\sec {\mathit{\theta }}_{\mathrm{0}}^{\star }+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{1}+\cos {\mathit{\xi }}^{\star })\sec {\mathit{\theta }}^{\star }\sec {\mathit{\theta }}_{\mathrm{0}}^{\star }.\end{array}$$

The term O in Eq. (A3) and the starred angles θ^⋆, ${\mathit{\theta }}_{\mathrm{0}}^{\star }$, and ξ^⋆ are calculated using the following set of equations:

$$\begin{array}{}\text{(A4)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}^{\star }=\arctan \left({\displaystyle \frac{b}{r}}\tan \mathit{\theta }\right),{\mathit{\theta }}_{\mathrm{0}}^{\star }=\arctan \left({\displaystyle \frac{b}{r}}\tan {\mathit{\theta }}_{\mathrm{0}}\right),\text{(A5)}& {\displaystyle }& {\displaystyle }\cos {\mathit{\xi }}^{\star }=\cos {\mathit{\theta }}^{\star }\cos {\mathit{\theta }}_{\mathrm{0}}^{\star }+\sin {\mathit{\theta }}^{\star }\sin {\mathit{\theta }}_{\mathrm{0}}^{\star }\cos (\mathit{\psi }-{\mathit{\psi }}^{\prime }),\text{(A6)}& {\displaystyle }& {\displaystyle }O={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}(t-\sin t\cos t)(\sec {\mathit{\theta }}^{\star }+\sec {\mathit{\theta }}_{\mathrm{0}}^{\star }),\text{(A7)}& {\displaystyle }& {\displaystyle }\cos t={\displaystyle \frac{h}{b}}{\displaystyle \frac{\sqrt{D^{\mathrm{2}}+{\left(\tan {\mathit{\theta }}^{\star }\tan {\mathit{\theta }}_{\mathrm{0}}^{\star }\sin (\mathit{\psi }-{\mathit{\psi }}_{\mathrm{0}})\right)}^{\mathrm{2}}}}{\sec {\mathit{\theta }}^{\star }+\sec {\mathit{\theta }}_{\mathrm{0}}^{\star }}},\text{(A8)}& {\displaystyle }& {\displaystyle }D=\sqrt{{\tan }^{\mathrm{2}}{\mathit{\theta }}^{\star }+{\tan }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{0}}^{\star }-\mathrm{2}\tan {\mathit{\theta }}^{\star }\tan {\mathit{\theta }}_{\mathrm{0}}^{\star }\cos (\mathit{\psi }-{\mathit{\psi }}^{\prime })}.\end{array}$$

The parameters $b/r$ and $h/b$ are the crown relative shape and the crown relative height, respectively. These were fixed to 1 and 2, respectively, following Strahler et al. (1999).

Appendix B: Dynamic gridding

For the version of the GOME-2 surface LER database described in Tilstra et al. (2017), v1.7, the spatial resolution was $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$. This is a relatively low spatial resolution, leading to artefacts, especially near coastal areas. In Fig. B1 this is illustrated by panels (a) and (b). Panel (a) presents the GOME-1 surface LER at 772 nm for March, for western Europe. Panel (b) presents the v1.7 GOME-2 surface LER. There are differences mostly because for GOME-2 the MODE-LER field was plotted but for GOME-1 the MIN-LER field was. The MODE-LER does a better job at detecting the snow-covered areas. Apart from these differences, both databases, having identical spatial resolutions, show similar difficulties near the coastline. The coastline seems to be pushed land-inwards. This is a direct consequence of the way the retrievals operate. By focusing on the smallest scene LER values collected in a grid cell, the observations over water are favoured over those over land.

Figure B1 Surface LER at 772 nm for the month of March from (a) the GOME-1 database and (b) the previous GOME-2 v1.7 database. Panels (c) and (e) help explain the new dynamic gridding procedure described in Appendix B. Panel (d) presents the GOME-2 database with coastline improvement, (f) does the same but with dynamic gridding also for mountain ranges, and (g) does the same but with the bilinear interpolation scheme applied. Panel (h) presents the Medium Resolution Imaging Spectrometer (MERIS) black-sky albedo as a reference for qualitative comparison. For the GOME-1 database the MIN-LER field was plotted and for the GOME-2 database the MODE-LER field was.

To remedy the artefacts, we first calculate the surface LER for three spatial resolutions: $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$, $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$, and $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$. The $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ field has the highest resolution and as such offers the best coastal representation. Unfortunately, the smaller grid cell size also results in less measurements per grid cell and in general leads to a lower quality mostly due to cloud contamination. In most cases, trading quality for a higher spatial resolution is not desirable. For coastal areas, however, the higher spatial resolution takes precedence.

We start with the $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ surface LER field and use it as a basis. We then perform a loop over all the $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid cells, and whenever the $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid cell does not contain a coastline, we overwrite the four associated $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ grid cells with the surface LER value of the overlapping $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid cell. Next, we loop over the $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid cells, and whenever the $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid cell does not contain a coastline, we fill the associated 16 $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ grid cells with the surface LER value of the overlapping $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid cell. The result is a database field that has an intrinsic resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ for most of the grid cells but an intrinsic resolution up to $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ near coastlines. This dynamic gridding procedure is illustrated in Fig. B1c for a part of the Portuguese coastline. Blue squares represent $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ grid cells, green squares represent $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid cells, and yellow squares represent $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid cells.

The coastline detection is performed using the Global Self-consistent, Hierarchical, High-resolution Geography (GSHHG) database (Wessel and Smith, 1996). The GSHHG database offers coastline information for the continents, islands, lakes, rivers, river-lakes, island-in-lakes, and even on the “pond-in-island-in-lake” level. We make use of the highest resolution available of the database, which is the “full resolution” version, available at https://www.soest.hawaii.edu/pwessel/gshhg/ (last access: 4 June 2021). For Antarctica we only consider the grounding coastline as a coastline. We do not take rivers and canals into account as these have negligible surface areas. We do take the so-called river-lakes into account. Islands with an area of less than 5000 km are not taken into account, nor coastlines from a “pond-in-island-in-lake”.

The improvement for coastal areas is demonstrated by Fig. B1d. Next, we focus our attention on the snow-covered mountain ranges in panel (d), which are captured poorly because of the low spatial resolution. For these regions, we manually assign rectangular parts of the grid as regions for which the intrinsic resolution should be fixed to $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$. In other words, in the dynamic gridding procedure these grid cells are protected such that they cannot be overwritten by the contents of the larger $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ and $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid cells. This procedure is illustrated in Fig. B1e. The regions containing mountain ranges are indicated in blue. The resulting surface LER field is shown in Fig. B1f. The Alps, Pyrenees, and Dinaric Alps are represented much better.

The approach described above leaves us with a database grid of $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ resolution that offers a higher resolution near coastlines and for snow-covered mountain ranges. However, for most of the regions over land and ocean the intrinsic resolution is $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$, and the surface LER grid is filled with redundant information in the form of 4×4 cell blocks in which 16 identical surface LER values are stored. This is not necessarily wrong, but it does complicate the interpolation that users need to apply to determine the surface reflectance for the measurements footprints they are dealing with.

To make things easier for the user the surface LER inside the 4×4 blocks is distributed over the 16 grid cells using standard bilinear interpolation over the 2D surface LER grid. Care is taken to only perform the interpolation inside and between grid cells that have an intrinsic resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ (i.e., the yellow grid cells in Fig. B1c and e). Other grid cells, such as the ones near the coastline or mountain ranges, are left untouched. After the bilinear interpolation a common additive correction factor is applied to the 16 grid cells in such a way that the average surface LER of these 16 grid cells is the same as before applying the bilinear interpolation. This step is needed because we do not want part of the reflectivity of the surface to disappear from the 4×4 blocks as a result of the bilinear interpolation. Next, we repeat the bilinear interpolation and apply the resulting additive correction factor to achieve a slightly higher level of smoothness, making the second 2D field a bit more convincing than the first one. The result is shown in Fig. B1g. It is important to stress that the above procedure is not an attempt to artificially increase the spatial resolution. It simplifies the interpolation that needs to be performed by the users.

Finally, in Fig. B1h we present the Medium Resolution Imaging Spectrometer (MERIS) black-sky albedo at 775 nm for comparison. The MERIS database has the same resolution as the GOME-2 surface LER database grid, so it can be used well for a qualitative comparison. Coastline and snow-covered mountain ranges compare quite well.