Surface Bidirectional Reflectance Distribution Function (BRDF) is a key property of the land surface describing its reflectance as a function of the sun-view geometry. The BRDF model was originally introduced by the MODIS and MISR science teams as standard products (Martonchik et al., 1998; Schaaf et al., 2002). The BRDF model has many uses: it defines surface albedo as a function of solar zenith angle (SZA) and aerosol optical depth (e.g., Lyapustin, 1999; Lucht et al., 2000) and provides albedo parameterization for use in climate models (e.g., Liang et al., 2005). It is used in Earth System Simulation Experiments (OSSE) to support mission formulation and sensor requirements and test processing algorithms. BRDF is widely used for normalization to a standard view geometry supporting detection of land surface change and fire-burnt scars (e.g., Roy et al., 2005; Giglio et al., 2018; Schläpfer et al., 2015), separating vegetation seasonality from variations in view geometry over tropics (e.g., Bi et al., 2015) etc. In the Multi-Angle Implementation of Atmospheric correction (MAIAC) algorithm, the 1 km spectral BRDF is used for accurate cloud and cloud shadow detection as well as for fast and un-ambiguous retrieval of snow fraction and snow grain size for detected snow (Lyapustin et al., 2018). BRDF-normalization became a standard component of on-orbit sensor calibration helping characterize long-term calibration trends and cross-calibrate different sensors over quasi-stable desert calibration sites (Lyapustin et al., 2014; Lyapustin et al., 2023).

The NASA MODIS and VIIRS data processing over land uses the reciprocal Ross-Thick Li-Sparse (RTLS) BRDF model introduced at the turn of the century (Lucht et al., 2000). This is a linear model using volumetric scattering and geometric-optics kernels which are functions of the sun-view geometry only. Linearity of this model allowed Lyapustin and Knyazikhin (2001) to develop an accurate semi-analytical expression for the top of atmosphere radiance as a function of three model parameters, and an accurate BRDF-coupled algorithm for atmospheric correction (Lyapustin et al., 2012; Lyapustin et al., 2018).

Over two decades of NASA MODIS observations have demonstrated a reliable performance of the RTLS model globally over snow-free land. Except near the principal plane, the RTLS model has a rather good accuracy over a snow as well although it falls behind the specialized two-parameter radiative transfer-based snow model with surface roughness (Lyapustin et al., 2010). Despite its general success, the RTLS model is not properly constrained at high zenith angles above 60°–70° (e.g., Gao et al., 2000). RTLS model works well for MODIS with its maximal scan angle of 55° which translates to the maximal view zenith angle (VZA) ∼62° when accounting for the Earth curvature. However, the VIIRS sensor, a follow-on to MODIS, has the maximal VZA∼72° at the edge of scan. The modern geostationary satellites (e.g., ABI on GOES 16–18; AHI on HIMAWARI 8–9) observe both the full range of solar zenith angle (SZA) variation, and high VZA near the disk edge. This stipulates the need to revisit the RTLS model and extend the range of valid zenith angles to ∼84° [μ = cos (VZA) ∼ 0.1].

RTLS is a semi-empirical linear model, represented as a sum of isotropic, geometric, and volume scattering components as originally introduced by Roujean et al. (1992): ρ μ 0 , μ , ϕ = k i + k v F v μ 0 , μ , ϕ + k g F g μ 0 , μ , ϕ . (1)

It uses predefined geometric functions (kernels), F _v , F _g , to describe different angular shapes. The kernels are independent of the land surface conditions. The bidirectional reflectance of a pixel is characterized by a combination of three kernel weights, K → = k i , k v , k g T .

The volumetric scattering (Ross-thick) kernel is expressed as: F v = μ ¯ + μ ¯ 0 − 1 π 2 − ξ cos ⁡ ξ + sin ⁡ ξ 1 + 1 + ξ ξ 0 − 1 − π 4 . (2)

Here, μ ¯ = μ and μ ¯ 0 = μ 0 are positively defined cosines of view and solar zenith angles and ξ is angle of scattering, cos ⁡ ξ = μ 0 μ + 1 − μ 0 2 1 − μ 2 cos ⁡ φ . The relative azimuth φ is defined such that the backscattering direction (hotspot) is at φ = 0 ° . Similar to Vermote et al. (2008), the V-kernel (Equation 2) is shown with the multiplicative hotspot factor with width of ξ 0 = 1.5 ° following Maignan et al. (2004), though we keep the original normalization of Lucht et al. (2000). We should mention that the standard RTLS kernel F _v does not include this hotspot term.

The geometric (Li-Sparse) kernel F _g is given by: F g = ∆ μ 0 , μ , ϕ + 1 + cos ⁡ ξ / 2 μ ¯ μ ¯ 0 , (3) ∆ μ 0 , μ , ϕ = 1 π t − sin ⁡ t ⁡ cos ⁡ t − π m , where m = 1 / μ ¯ + 1 / μ ¯ 0 is an airmass factor, and cos ⁡ t = 2 / m t g 2 θ 0 + t g 2 θ − 2 t g θ 0 t g θ ⁡ cos ⁡ φ + t g θ 0 t g θ ⁡ sin ⁡ φ 2 .

The geometric term represents reflectance reduction by shadows from the elevated ellipsoidal protrusions (e.g., trees), and thus the general hot spot effect (Li and Strahler, 1992). The volumetric kernel represents light scattering in vegetation media with randomly oriented leaves. It is a linearized version developed by Roujean et al. (1992) based on the single scattering model in vegetation Ross (1981).

The RTLS kernel shapes are illustrated by black lines in Figure 1 in the principal plane (0°–180°), cross plane, and diagonal plane (45°–135°) which is close to a typical view geometry of MODIS.

Both F _v at SZA = 84.3° and F _g at all shown solar angles rapidly increase in amplitude past 60°–70° and are truncated in the plots before reaching maximal VZA = 84.3°. This instability of kernels is generated by μ ¯ , μ ¯ 0 or their sum in denominator and by the air mass factor m, specifically μ ¯ + μ ¯ 0 − 1 in F _v and in terms m and 1 / 2 μ ¯ μ ¯ 0 in F _g . Limitations of the RTLS model at high zenith angles were discussed in the relevant publications (e.g., Li and Strahler, 1986; Wanner et al., 1995): they include the lack of multiple shadowing, assumption of a black solar and view (obscuration) shadows, the equivalency of the solar and obscuration shadows, the neglect of the diffuse light, of effect of the tree/canopy height variations and others. Our goal is not to develop a much more complex BRDF model by incorporating these effects, but rather suggest a simple normalization or scaling solution which improves model performance at high zenith angles. For practical purposes, we consider the range of angles up to ∼84° (μ ∼ 0.1). The idea is to slow down a rapid increase of μ − 1 at VZA > 60°. For instance, replacing μ with its square root at μ = 0.1 gives 1 / μ ∼3.16 instead of 1 / μ = 10 . Thus, the barred symbols are defined as linearly scaled between μ and μ below μ = 0.5 : μ ¯ = μ , μ ≥ 0.5 , μ ¯ = w μ + 1 − w μ , w = μ / 0.5 , μ < 0.5 . (4)

In summary, the scaled sRTLS model, given by Equations 1–4, differs from the standard RTLS by adding the hotspot factor in F _v and using scaled values μ ¯ , μ ¯ 0 at zenith angles above 60°. Scaled kernels, illustrated in red in Figure 1, show constrained behaviour at high zenith angles. The difference increases with SZA and is highest at SZA = 84.3°.

MAIAC algorithm has long been adapted to processing of HIMAWARI-8,9 AHI and GOES 16–18 ABI observations (e.g., Li et al., 2019; She et al., 2019) with some further research developments (Wang et al., 2022). MAIAC sits at the core of the GeoNEX (Geostationary NASA Earth eXchange) project, a collaborative effort led by scientists from NASA, NOAA, and many other research institutes around the world (Nemani et al., 2020). The goal of GeoNEX is to create a suite of globally unified and consistent science products from the constellation of latest operational geostationary satellite sensors (Chen et al., 2021; Wang et al., 2020; Hashimoto et al., 2021; Wang et al., 2022; Zhang et al., 2022; Li et al., 2023; Shen et al., 2023; Gao et al., 2024). An atmospheric correction algorithm that rigorously accounts for varying sun-sensor geometries and surface BRDF is thus essential to the success of all the downstream products and the overarching goal of GeoNEX.

The original MAIAC-Geo includes standard MAIAC algorithm components (Lyapustin et al., 2018), namely, cloud, cloud shadow and snow detection, aerosol retrieval and atmospheric correction producing spectral bidirectional reflectance factors (BRF), commonly called surface reflectance, spectral albedo and daily BRDF model parameters. Similarly to MAIAC MODIS, we also provide snow fraction and snow grain size for detected snow, and smoke plume injection height for detected wildfires (Lyapustin et al., 2020).

For each pixel, the MAIAC-Geo atmospheric correction uses the BRDF-coupled radiative transfer model (Lyapustin et al., 2012; Lyapustin et al., 2018) and the BRDF model from the previous day to derive a BRF for each ABI observation acquired every 10 min. At the end of the day, it combines all daily BRFs and derives an updated BRDF model that provides the best fit to the multi-angle BRF dataset. At the onset, when BRDF is not initialized, MAIAC uses a Lambertian surface reflectance model and computes Lambert Equivalent Reflectance (LER). Once an initial approximation of the BRDF model is derived from LERs, the full BRF/BRDF processing rapidly converges to its full accuracy in 1–2 days. The full BRDF inversion, producing parameters k i , k v , k g } is performed on days with low-to-moderate cloud coverage, when observations cover the full range of SZA and scattering angles. On partly cloudy days, with at least 10 cloud-free observations and representative angular coverage of low, medium and high SZA ranges, we adjust the BRDF magnitude by scaling it to the retrieved BRFs assuming stable BRDF shape from previous retrievals. This approach is similar to the MODIS BRDF\albedo backup algorithm, flagged as lower quality magnitude inversion (Schaaf et al., 2002), and is also used in the MAIAC MODIS/VIIRS algorithms (Lyapustin et al., 2018). In both inversion and scaling cases we calculate the standard deviation (σ) representing daily goodness of fit.

This work uses GOES-16 ABI data over the East Coast United States in 2020 and 2023 to illustrate performance of both the standard RTLS and proposed sRTLS models. The annual variation of the view geometry (scattering angle) is illustrated in Figure 2. For most pixels, excluding those near the edge of disk at extreme view angles, the scattering angle approaches the hotspot (exact backscattering) twice a year according to seasonal change in the local solar elevation (Li et al., 2021). The total range of daily variation of scattering angle changes little between second decades of March and October but decreases towards the turn of the year. We should mention that MAIAC follows the convention used in the radiative transfer theory where the backscattering corresponds to the relative azimuth and scattering angle of 180° contrary to the definition used in the RTLS model.

Figures 3, 4 show the results of processing for two pixels representing an urban area of Washington, DC (Pix1), and a forested area about 100 km north-east (Pix2), viewed approximately at VZA ∼ 45°. In both cases, we show the results typical for the hotspot view geometry (a) and for the regular geometry (b). The dots show the BRF retrieved using RTLS (black) and sRTLS (red) BRDF models, and the LER (blue). The empty circles connected by solid line show the best fit BRDF models. The first Figure 3A shows 5 consecutive days 64–68 of 2023 of mostly cloud-free observations. During this time period, SZA_max ≈ VZA which captures the hotspot with maximal daily scattering angles between 177.5° and 179.1°. The retrievals are shown for four channels, Blue (0.47 μm), Red (0.64 μm), NIR (0.86 μm) and shortwave IR (SWIR, 2.24 μm). Both sRTLS and RTLS models provide a very similar match to the retrieved BRFs. The main differences appear near the hotspot and at high SZAs. Due to the Maignan hotspot factor, the sRTLS model fits the hotspot direction well: the BRF-BRDF difference near the hot spot does not exceed 0.01 in all four channels. The RTLS model significantly underestimates the hotspot reflectance by 0.025 (Blue), 0.04 (Red), 0.055 (NIR), 0.052 (SWIR). At high SZAs, the sRTLS model gives higher values of BRF and by design, higher BRDF. In the legend, the BRF obtained with sRTLS model is called sBRF. The difference sBRF-BRF as well as sBRF-LER is negligible in the SWIR but it increases towards the Blue with increase of the diffuse/direct radiation ratio due to higher atmospheric optical depth and Rayleigh-aerosol scattering of light. For instance, at SZA ∼ 80° the difference reaches sBRF-BRF ∼ 0.02 in the blue channel. Such difference may persist throughout the year (e.g., Figure 3B).

The difference in BRDF shape between the two models at high SZA is strongly manifested for the forested pixel (Pix2). Figure 4A shows the reflectance values for days 284–286 coinciding with the fall hotspot period (max γ = 178.7°). Similarly to the urban pixel (A), the RTLS model underestimates the reflectance from 0.012 (Blue) to ∼0.05 (NIR), while the sRTLS model provides the good match to the BRF. The RTLS BRDF displays a characteristic “shelf” at SZA ∼ 58°–65° with rapid decrease at higher angles. The difference in data fit is best visible in the SWIR channel. In this case, sBRF and BRF are nearly identical, and the sRTLS model provides a significantly better match.

Figures 3, 4 also shows the LER by blue dots. The Lambertian surface model is widely used for atmospheric correction, e.g., by the MODIS and VIIRS standard AC algorithm (Vermote and Kotchenova, 2008). This approximate model simplifies the radiative transfer problem and allows for faster atmospheric correction using significantly smaller look-up tables. The LER has a good accuracy at longer NIR-SWIR wavelengths, when Rayleigh and aerosol optical depths are low, but it has a systematic bias at shorter wavelengths (e.g., Wang et al., 2010; Lyapustin et al., 2021) which is also apparent in Figures 3, 4. The maximal LER bias ∼22% is near the hotspot signature.

Figure 5 shows histograms of daily standard deviations (σ) from all pixels over an area of 250 × 250 km² centred at the NASA Goddard Space Flight Center, during the full year of 2020. The results are shown for four channels, with Blue (0.47 μm) and Red (0.64 μm) represented by the respective colors, and NIR (0.86 μm) and SWIR (2.24 μm) shown in brown and black, respectively. Solid and dashed lines represent the sRTLS and the standard RTLS models, respectively. One can see that sRTLS model gives a better fit in all bands shifting histograms of standard deviation towards smaller values. The use of sRTLS model slightly improves the mean standard deviation, computed as the normalized first moment of the distribution, by (2-4) 10^–4, in the Blue-SWIR bands except Red (10^–5). While the mean reduction of σ is low, being distributed over up to 60–100 observations per day, the main impact on both BRFs and BRDF fit is at high zenith angles and near the hotspot and can be significant as we demonstrated earlier.

Developed in late 1990s, RTLS has been a highly successful model describing BRDF of the land surface for climate models and many applications listed in the Introduction. The model works well for zenith angles up to ∼60°–70° but its accuracy deteriorates at higher zenith angles and near the hotspot. This limits its use in processing of geostationary observations and in applications using high zenith angles. In this work, we described a simple practical solution that allows to extend the range of the sRTLS model to higher zenith angles. sRTLS improves the accuracy both near the hotspot directions and at zenith angles above ∼70° which is achieved by adding the multiplicative hotspot factor of Maignan et al. (2004) and by scaling the cosines of zenith angles in the RTLS model at ZA > 60°, respectively. We tested the model using GOES-16 ABI data processing for zenith angles up to ∼80° which is currently set as a limit in MAIAC-Geo look-up tables. Further extension to higher zenith angles would require consideration of the Earth sphericity (e.g., Korkin et al., 2022) in MAIAC look-up tables which is currently under development. It should be mentioned that the introduced modification of the BRDF model is expected to have a very minor effect on the surface albedo except at high SZAs where the difference may become noticeable. The analysis related to surface albedo will be published elsewhere.

In MAIAC, knowledge of surface BRDF is used for detection of clouds and cloud shadows. At high solar zenith angles, both problems increase in complexity as the contrast from clouds diminishes and the length of shadows grows. Thus, the improvement in the BRDF model performance at high zenith angles becomes critical for the total processing accuracy of geostationary data, including both aerosol retrieval and atmospheric correction. In this regard, the particular benefits of the sRTLS model are expected at high latitudes which experience the most rapid climate change effects. Through the improved aerosol data quality in the morning and pre-dawn hours, this development is also expected to benefit the air quality monitoring from geostationary satellites, in particular the analysis of diurnal cycle of air pollution related to traffic commute patterns, industrial activities and other sources.

The sRTLS model has been implemented in MAIAC MODIS Collection 6.1 algorithm, in version 2 MAIAC VIIRS algorithm, and in MAIAC-Geo algorithm for AHI HIMAWARI and GOES-16–18 ABI data. The spectral BRF, sRTLS model parameters and albedo from ABI and AHI are reported in the NASA GeoNEX Level 2 land products, which are publicly accessible from the project data portal (https://data.nas.nasa.gov/geonex/data.php). In addition, a set of open-source software tools (in R and Python) has been developed to subset, analyse, and visualize these products at specified locations and time ranges, which are also freely accessible online (https://github.com/nasa/GeoNEXTools).