The brightness test aims to detect optically thick bright clouds that have a high albedo in the visible spectrum. A pixel is masked as cloud if the measured reflectance (R^m) exceeds the theoretical value at maximal AOD = 6 of the MAIAC look-up table (LUT) at a given view geometry (R _max) with a certain threshold: R m > R max + Thresh (1)

The threshold is 0.1 over bright Sahara region and 0.05 otherwise. The brightness test uses the EPIC blue channel (443 nm), where the surface is generally dark, and the reduction of TOA reflectance by absorbing aerosols (smoke/dust) and the difference in reflectance with non-absorbing clouds is maximal compared to longer wavelengths.

In general, clouds exhibit a larger spatial variance than aerosol and the ocean surface. The spatial variance test computed using 3 × 3 pixel window has been a standard technique for cloud detection over the ocean at moderate resolution ∼1 km (e.g., Martins et al., 2002). Based on simulated EPIC observations from 1 km MODIS data, we selected the spatial variance threshold of 0.005 which achieves a reasonable balance between cloud filtering and fraction of clear pixels.

Over land surfaces, using a global threshold is problematic due to spatial variability of the land surface reflectance, in particular over bright deserts, in the urban regions and over agricultural areas. Working with gridded data, MAIAC keeps memory of the 3 × 3 standard deviation (σ) for each 10 km grid cell derived in cloud-free and low aerosol conditions. Similar to MAIAC MODIS, σ is computed for the red and blue bands and updated on cloud-free days from observations closest to nadir, when the observation footprint is minimal and spatial variance from surface is maximal. The implementation for EPIC follows test (C.4) in Lyapustin et al. (2018).

Detection of optically thin high clouds relies on EPIC measurements in the oxygen A-band. While this signal is low in cloud-free conditions due to absorption by molecular oxygen, presence of high clouds creates a relatively strong signal. Detection of high clouds employs a reflectance ratio measurement of oxygen A-band to the window channel (780 nm) divided by a theoretical reflectance ratio for surface elevation (Z) (Zhou et al., 2020): R 764 m e a s / R 780 m e a s R 764 t h e o ( Z ) / R 780 t h e o ( Z ) = exp ( − mτ O 2 A ) exp ( − m τ O 2 A , Z ) (2) where τ_O2A and τ_O2A,Z are optical depth values due to O₂ absorption from the measured ratio and from the theoretical LUT, respectively, and m = 1/μ ₀ + 1/μ is an atmospheric airmass factor depending on cosines of solar (μ ₀), and view (μ) zenith angles. A simple threshold-based approach is then used for high cloud detection: the pixel is considered cloudy if τ O 2 A , Z − τ O 2 A > 0.056.

Similarly to MAIAC MODIS, the above tests only serve for an initial cloud screening. The cloud mask is significantly enhanced following aerosol retrievals by limiting the small-scale spatial variability of AOD, and during the atmospheric correction through comparison of spectral reflectance with the predicted values based on the BRDF model.

Following MAIAC MODIS Collection 6 algorithm, we are using eight prescribed regional aerosol models to represent variability of aerosol properties over global land. Geographic distribution and model parameters are provided in Lyapustin et al. (2018, Figure 4 and Table 1).

One known issue in MAIAC C6 was underestimation of AOD for the biomass burning aerosol at high AOD (e.g., Lyapustin et al., 2018; Schutgens et al., 2020; Sogacheva et al., 2020). As a remedy, in MODIS MAIAC C6.1 we adjusted the model parameters at AOD>0.6 based on the regional climatology analysis of the AERONET (Holben et al., 1998; Giles et al., 2019; Sinyuk et al., 2020) record. Similarly, to correct the known low bias of the mineral dust AOD over Western Sahara, we introduced a new corresponding region with the more absorbing dust model. These amendments are used in the v2 MAIAC EPIC and will be described in detail elsewhere.

MAIAC processing uses the ancillary NCEP ozone and column water vapor information. The over ocean processing also uses the NCEP wind speed.

Retrieval of SRC is a central component of MAIAC: it provides separation of the surface and atmospheric signals in the TOA measurements, and is required for aerosol retrievals. Because EPIC lacks the 2130 nm channel used in MAIAC MODIS, we define SRC as the ratio of the surface reflectance in Blue to Red (SRC = ρ_L,Blu/ρ_L,Red) bands. Reflectance ρ_L,λ is a Lambertian reflectance resulting from the Rayleigh atmospheric correction with low background aerosol. SRC is obtained as a minimum value over the 2-month period of time. Following Lyapustin et al. (2018), we are using two independent lines of update, shifted by 1 month. Thus, the SRC is dynamically updated at least every month or more frequently if a new minimum value is found. SRC is characterized in 4 bins in cosine of the solar zenith angle 1–0.9, 0.9–0.7, 0.7–0.45, and 0.45–0.2. The SRC for the morning and the afternoon observations is separate because of the change in geometry at a near-constant scattering angle, e.g., depending on the part of the orbit, the view zenith angle (VZA) may be higher than the solar zenith angle (SZA) in the morning but symmetrically lower in the afternoon and vice versa.

The AOD is obtained by matching the observed and theoretical TOA reflectance at 443 nm based on the look-up table. The surface reflectance ρ_B is evaluated from the atmospherically corrected ρ_L,Red (AOD) at 680nm, ρ_L,Blu = SRC*ρ_L,Red. This AOD is derived using the corresponding regional background aerosol models. When derived AOD is high (>0.6) and absorbing smoke or dust is detected, the v2 MAIAC runs a separate inversion of UV-vis observations providing AOD and spectral aerosol absorption (Lyapustin et al., 2021a, this issue).

At high altitudes (over 3.5 km, e.g. Tibetan plateau) where AOD is generally very low and MAIAC AOD retrievals do not have sufficient accuracy, we assume a fixed climatology AOD _min = 0.02 for the atmospheric correction.

To assess accuracy of AOD retrieval, we conducted AERONET validation for 2015–2020 using level 2.0 AERONET version 3 data (Giles et al., 2019). The comparison uses EPIC AOD at 443 nm averaged over the 5 × 5 pixels window (50 km) with AERONET data selected within 30 min from the satellite observation. The EPIC data were filtered according to the Sun and view zenith angles less than ∼63° and at least 50% coverage in the spatial window.

The validation results are presented in Figure 2. The site-level global statistics shows the correlation coefficient (R), root mean square error (RMSE) and the mean bias error (MBE, MAIAC–AERONET). MAIAC shows a good retrieval accuracy, with R ≥ 0.7–0.8 and low RMSE and MBE, over vegetated parts of the world including North and South America, north-central Eurasia and Oceania. AOD is generally overestimated over bright surfaces such as western United States and Australia. Such bias is typical for aerosol products based on a single-view satellite measurement, and it is exacerbated for EPIC due to unfavorable view geometry near the backscattering direction. The underestimation of AOD is obvious in regions of strong biomass burning, including central Africa, Indo-Gangetic plain and south Asia. Similar to the bias, the RMSE is generally low globally with the exception of the major dust and biomass burning aerosol source regions which also have a much higher annual average AOD. The correlation is generally high in regions with higher magnitude and variability of AOD. On the contrary, the low correlation is observed over regions with low aerosol loading and low natural variability, such as Australia or south of the South American continent.

The two scatterplots on the right summarize the global validation analysis. The left one 1) shows all AERONET sites. On the right scatterplot (b), we excluded 11 bright surface sites over the Western United States (Bakersfield; Goldstone; KeyBiscayne; Neon_ONAQ; Railroad Valley; Sandila_NM_PSEL; TableMountain_CA; Tucson; UACJ_UNAM_ORS; White_Sands_HELSTF; Yuma) and one site over Australia (Birdsville). These sites are located in arid regions with low AOD and low AOD variability where MAIAC EPIC overestimates AOD and shows low R-values. Due to low cloudiness, these sites also contribute a disproportionate ∼9% of the total matching points. Considering the right scatterplot as a baseline, MAIAC shows an overall correlation of 0.77 with RMSE ∼ 0.159 and a mean bias of 0.046. Over 66.72% of MAIAC EPIC AOD agree with AERONET within the expected error (EE) of ±(0.1 + 0.1AOD). The v2 shows an improvement over v1 which had a global statistics of R = 0.69, RMSE = 0.17, MBE = 0.03 (unpublished).

The v1 MAIAC EPIC algorithm used a Lambertian model for the atmospheric correction, where the surface reflectance is derived from the following approximation to the TOA reflectance: R T O A ( μ 0 , μ , φ ) ≈ R A ( μ 0 , μ , φ ) + ρ T d ( μ 0 ) T u ( μ ) / ( 1 − s ρ ) , ρ ≡ ρ ( μ 0 , μ , φ ) (4)

Eq. 4 only requires the knowledge of atmospheric (path) reflectance (R^A), upward (T^u) and downward (T^d) atmospheric transmittance as functions of the cosines of Sun (μ ₀) and view (μ) zenith angles and relative azimuth (φ), and spherical albedo of atmosphere (s). In the EPIC view angles near the hotspot where surface is brighter than in the other directions, the Lambertian approximation underestimates the surface reflectance (e.g., Wang et al., 2010). An analysis of Lambertian biases was recently given by Lyapustin et al. (2021b) based on a comparison between the two MODIS surface reflectance products, the standard surface reflectance (SR) MOD09 based on Lambertian assumption, and the bidirectional reflectance factors (BRF) of algorithm MAIAC (MCD19A1).

The MAIAC atmospheric correction uses a rigorous expression for the TOA reflectance. Taking advantage of linearity of the RTLS function, it represents the TOA reflectance explicitly using weights of the RTLS model: R T O A ( μ 0 , μ , φ ) = R A ( μ 0 , μ , φ ) + k L F L ( μ 0 , μ ) + k V F V ( μ 0 , μ , φ ) + k G F G ( μ 0 , μ , φ ) + R n 1 ( μ 0 , μ ) (5)

Here, F-functions are integrals of the atmospheric path radiance incident on surface and atmospheric Green’s function (Lyapustin and Knyazikhin, 2001) with respective kernels of the RTLS model. R^nl is a weakly non-linear function of the surface reflectance, describing multiple light scattering between the surface and the atmosphere. The F-functions and R^nl are computed analytically using eight primary functions which are stored in the MAIAC LUT (e.g., Lyapustin et al., 2011a). For the purpose of atmospheric correction, let us re-write Equation 5 as follows: R T O A ( μ 0 , μ , φ ) = R A ( μ 0 , μ , φ ) + c R S u r f ( μ 0 , μ , φ ) (6) where R^Surf is a surface-reflected term combining the last four terms of Eq. 5, and c ≡ c _λ is a spectrally dependent scaling factor. The R^Surf is computed using the BRDF model parameters stored in MAIAC memory for each grid cell (e.g., 1 km for MODIS and 10 km for EPIC). Then, the BRF is given by a scaled value: r λ = c λ R T L S λ ( μ 0 , μ , φ ) (7) where RTLS_λ is the BRDF model value for a given geometry. Because R^Surf is a nonlinear function of the surface reflectance, solving eqs. (6) and (7) takes 2 iterations (see Lyapustin et al., 2018, p. 5753).

Implementing rigorous atmospheric correction given by Eqs. (6) and (7). requires knowledge of the entire BRDF shape to correctly represent surface reflection of the direct Sun beam and of the diffuse (sky) light. In v2 MAIAC EPIC, we used the global Collection 6 MAIAC MODIS 1 km BRDF product MCD19A3 to develop an ancillary BRDF dataset for the EPIC 10-km grid from the closest MODIS channels. The AC based on Eqs. (6) and (7). uses scaling approach and only requires knowledge of the relative BRDF shape rather than the absolute reflectance. For this reason, the wavelength difference between the EPIC and MODIS channels is not important as long as the land surface reflectance among the paired channels, and thus the BRDF shape, remain similar. The ancillary 10 km BRDF for EPIC was created for every month starting in 2015. Figure 3 illustrates the global RGB BRDF for nadir view and Sun at 45° for January and July of 2016.

Figure 4 gives an example of atmospheric correction based on scaling (S) in the Blue and Red bands in columns 3 and 5, respectively. The result is shown as an excess of anisotropic over the Lambertian reflectance (BRF-ρ_L)/ρ_L (%). In agreement with theory, the difference is lowest when the Sun is near zenith, and it grows with the atmospheric airmass factor. It also increases with the total atmospheric optical thickness from NIR to Blue, for instance from ∼3% at 780nm, ∼10% at 680nm to 60–80% at 443 nm at high Sun/view zenith angles (μ ₀,μ ∼ 0.2). While this pattern agrees with theoretical expectations in general, the strong increase of BRF at high zenith angles, in particular at 443nm, does not seem realistic.

The scaling approach Eqs. (6) and (7) assumes that the BRDF model gives a good description of the surface reflectance at the angles of satellite observations. In this case, direct and diffuse surface-reflected signals at the TOA can be scaled using the same multiplier c. This approach works for MAIAC MODIS where the BRDF represents the view geometry sampled by MODIS. The ancillary MODIS BRDF was derived for the range of SZA observed near the local noon around 10:30 am (Terra) and 1:30 pm (Aqua) equatorial crossing time. Thus, it can be considered representative of the EPIC’s view geometry within about ±2 h of the local noon. Outside of that range, at higher SZA both earlier in the morning and later in the afternoon, the MODIS BRDF can still be used to compute the reflection of the diffuse sky irradiance assuming BRDF reciprocity, at least for the range of SZA agreeing with the VZA range of MODIS, ∼0–62° accounting for the Earth’s curvature. However, it cannot represent correctly the direct TOA reflectance ρ(μ ₀ ,μ,φ)exp(−mτ), where τ is an atmospheric optical depth, as MODIS does not make measurements at higher SZAs near the hotspot.

In this case, we can single out the direct reflectance in Eq. 6: R T O A ( μ 0 , μ , φ ) = R A ( μ 0 , μ , φ ) + ρ ( μ 0 , μ , φ ) e − m τ + g λ R D i f ( μ 0 , μ , φ ) (8) where the diffuse component of the surface-reflected signal at TOA is: R D i f ( μ 0 , μ , φ ) = R S u r f ( μ 0 , μ , φ ) − R T L S ( μ 0 , μ , φ ) e − m τ (9)

Above, ρ is the true surface BRF, and R ^Dif and RTLS are computed with the ancillary MODIS BRDF. The g_λ is the spectral adjustment factor designed to account for the surface reflectance difference from the spectral shift between the paired EPIC - MODIS channels, e.g., for the Blue (443/465.5 nm), Green (555/553.5 nm), Red (680/644.9 nm) and NIR (779.5/855.6 nm) EPIC/MODIS center wavelengths, respectively. For each 10 km pixel, we compute the spectral adjustment factor g_λ using scaling atmospheric correction (Eqs. (6) and (7) near the local noon (within Δμ₀ of ±0.1) where MODIS BRDF should be valid for both direct and diffuse terms, as we discussed above. In this case, g_λ is equivalent to the scaling factor c_λ. The g_λ is computed daily for each cloud-free 10 km grid cell at low AOD and is stored in MAIAC memory until updated with the next retrieval. Such approach allows us to evaluate the diffuse reflected term using the ancillary MAIAC MODIS BRDF, and compute BRF (ρ) from the direct reflected term in Eq. 8.

The described “direct term” (D) algorithm is more generic than the scaling approach. The resulting BRF_D (columns 4 and 6 of Figure 3 for the Blue and Red bands, respectively) can be compared to the scaling BRF_S in Figure 4. BRF_D shows a more constrained increase over the Lambertian reflectance up to SZA∼70° which grows only to ∼35% in the Blue band instead of ∼60–80% for the scaling approach. Importantly, the anisotropic enhancement (of the Lambertian reflectance) remains nearly constant in the range of SZA∼0–60°, though the uncertainty Δρ increases at high EPIC SZA/VZA. The uncertainty is twofold: it is related to both the aerosol retrieval uncertainty (Δτ) and to the uncertainty in the diffuse signal R ^Dif introduced by the MODIS BRDF (ΔRTLS) which was defined for a relatively narrow range of SZA values of MODIS Terra and Aqua at a fixed overpass time, and the limited range of VZA≤62°. Moreover, while working well in the range of SZA/VZA∼60°, the RTLS BRDF model has the problem of unconstrained growth of both geometric-optics and volumetric kernels at higher SZA/VZA, proportionally to a combination of terms 1/μ ₀, 1/μ (e.g., Gao et al., 2000). Finally, a significant uncertainty is related to the increase of the EPIC’s footprint with VZA faster than 1/μ while the ancillary BRDF generated from 1 km MODIS still closely represents the 10 km grid box. For these reasons, the atmospheric correction problem (Eqs. (6) and (7) at high zenith angles becomes ill-posed and poorly constrained, with an exponential propagation of uncertainties: Δ ρ ( μ 0 , μ , φ ) ∼ − e m τ { Δ R A ( Δ τ ) + g λ R D i f ( Δ τ , Δ R T L S ) } (10)

The above analysis showed limitations of both scaling and the “direct term” atmospheric correction methods, in particular at high SZA/VZA. The main limitations stem from the limited angular sampling of EPIC prohibiting deriving the self-consistent BRDF model in the full hemisphere of angles of incidence and reflection, and from the growing uncertainties of the ancillary MODIS BRDF model at high SZA/VZA in application to EPIC. Both Lambertian and scaling algorithms reproduce well the spatial pattern and the RGB color of the EPIC TOA images while, respectively, underestimating and overestimating the true BRF, especially at high zenith angles. The “direct term” algorithm shows rapidly growing uncertainties at high zenith angles. As it depends on the absolute ancillary BRDF model, this approach is also prone to spatial and spectral distortions in the resulting RGB BRF images. At the same time, this algorithm provides a realistic more constrained BRF increase over the Lambertian value, and a near-constant u = BRF/ρ_L ratio in the wide range of zenith angles up to ∼50–60°. This ratio fully agrees with the BRF/ρ_L ratio of the scaling method evaluated for the observations near the local noon.

Given these findings, the MAIAC EPIC v2 AC approach in RGB and NIR bands is implemented as follows:

- Compute the Lambertian reflectance (ρ_L,λ) from Eq. 4;

The anisotropic conversion factor near the local noon gives the increase over the Lambertian reflectance from 1–2% at 780 nm to 8–15% at 443 nm. The uncertainty of the reported BRF is low for the EPIC observations near the local noon, and it is expected to significantly increase at zenith angles above ∼60°. The selected empirical AC approach is fast and does not create spectral distortions, but it probably underestimates BRF at higher zenith angles. In near future, we plan to further explore both “scaling” and the “direct term” AC algorithms using BRDF model from the geostationary satellites which provide the full range of the solar zenith angle variations.

The MAIAC v2 algorithm reports both ρ_L and BRF reflectance values in the RGB and NIR channels. In the UV, where uncertainties are the largest, only the Lambertian reflectance is reported.

Figure 5 gives an example of atmospheric correction for the 1-month period of June 2–July 2, 2020, for a single 10 km bright surface grid cell in Arizona, United States. The derived surface reflectance displays a well-reproducible daily pattern in the visible–near infrared with surface reflectance increasing with SZA. The strongest growth is observed in the NIR in agreement with (Marshak, 2021). The pattern becomes less certain at UV wavelengths: the SR points tend to cluster within about ±0.01 of the mid-day value with occasional high and low outliers at both 388 and 340 nm wavelengths. On average, the AC produces a correct pattern with ρ₃₄₀<ρ₃₈₈ at μ,μ₀ above ≈0.6 (53°), and unstable result at higher SZA/VZA. This result is rather systematic and holds over both bright and dark vegetated surfaces. As the uncertainty in the UV SR rapidly grows with the airmass factor, the most accurate SR values at 340 and 388 nm are reported near the local noon.

Following computation of BRF, MAIAC proceeds with calculation of the RTLS parameters (K-coefficients) using the multi-angle dataset accumulated in the MAIAC Queue memory for each grid cell (up to 40 observations). This retrieval is performed for the four visible and near-IR bands only. After inversion, the new BRDF goes through several tests verifying “correctness” of its shape, and its consistency with the previous solution (Lyapustin et al., 2012). Figure 6 shows the daily BRF pattern (dots connected by solid lines) in the Blue-NIR EPIC bands for the pixel displayed in Figure 5 on two different days in June of 2020. The best-fit BRDF model is shown by the dashed lines. While the BRDF model error can reach several absolute percent of reflectance at high zenith angles, the typical rmse of the fit is low, within ∼0.001–0.002.

It should be mentioned that the derived BRDF model describes only the range of the EPIC observations within 4–12° from the backscattering direction, and is not representative of the general BRDF shape. A failure of the DSCOVR gyroscopes in late June of 2019 led to EPIC being placed in a safehold mode till March of 2020. This period was required to find the engineering solution for the satellite navigation using startrackers and update the geolocation algorithm. After resuming operation, the DSCOVR orbit became less constrained and allows the range of angles ∼2–12° from the exact backscattering since March 2020.