1.1.

Background

Following the recommendation of the Committee on Earth Observation Satellites (CEOS) Working Group on Calibration and Validation (WGCV), the Libya-4 PICS, which is one of the driest, least cloudy, and highly reflective desert PICS,¹² has been used extensively for post-launch radiometer drift monitoring and sensor pair radiometric scaling of high, medium, and low-resolution satellite imagers, including Landsat, Sentinel-MSI/OLCI/SLSTR, AVHRR, ATSR, MERSI, MODIS, VIIRS, DSCOVR-EPIC, and Meteosat GEO imagers. Libya-4 is located in an extensive sand dune region nearly void of vegetation, positioned near the northeast Libya and Egypt border, with an elevation of 118 m and monthly mean precipitation of $\sim 1\mathrm{mm}$. Staylor¹³ discovered the potential of the nearby Saharan desert PICS (17°N to 27°N and 3°W to 27°E) for monitoring the CERES precursor, Earth radiation budget experiment (ERBE) broadband SW channel calibration stability. Later, Staylor¹⁴ was able to determine the degradation of the AVHRR sensors utilizing the Libyan desert (20°N to 30°N and 20°E to 30°E) based on the ERBE calibration approach. Rao et al.¹⁵ utilized the Libyan desert (21° to 23°N and 28°E to 29°E) to monitor the degradation of the NOAA-7, 9, and 11 AVHRR sensors. Cosnefroy et al.¹⁶ established the Saharan and Arabian PICS locations, nomenclature, and initial stability assessment, including the Libya-4 site (28.55°N and 23.39°E), that were suitable for assessing the visible sensor calibration over multiple solar and viewing conditions. Over the years, follow-on studies have confirmed that Libya-4 is indeed one of the most stable and spatially uniform PICS suitable for sensor stability validation.^17–19 Due to the stability of the Libya-4 PICS, it was utilized to tie the 40-year AVHRR record to the Aqua-MODIS^20–22 and SeaWiFS²³ calibration references.

Once the temporal stability of Libya-4 was established, the observed viewing and solar angular dependencies of the surface BRDF were characterized. The Libya-4 sand dunes have an average height of 60 m and tend to reflect more in the backscatter than in the forward scatter direction, especially for large solar angles.²⁴ Bacour et al.¹⁹ stated that the Libya-4 multi-angle PARASOL observations reveal a surface anisotropy of less than 10%. Wu et al.²⁵ found that the Libya-4 BRDF is nearly Lambertian for viewing zenith angles (VZA) less than 25 deg and for solar zenith angles (SZA) less than 50 deg by comparing Lambertian and semi-empirical BRDF models. Bhatt et al.²⁶ examined Meteosat-9 imager daily clear-sky TOA radiances under off-nadir (VZA of $\sim 42\mathrm{deg}$) viewing conditions over Libya-4 and reported that the average daily inter-annual variability for visible wavelengths was below 1%. Given the Libya-4 minimal interannual variability, the CERES IGCG radiometrically scaled successive Meteosat geostationary imagers operating at 0 deg longitude to a common calibration reference.²⁷ Chang et al.²⁸ found that the Terra and Aqua MODIS derived semi-empirical and empirical BRDF models were similar. Mishra et al.²⁹ derived a Terra-MODIS TOA reflectance model with respect to SZA and remarked that higher order Ross–Li kernel functions were no better than the simple linear model for near nadir measurements.

Because the BRDF effects are small over desert sites, many desert calibration studies have simply correlated sensor reflectances with SZA. Wu and Zhong³⁰ utilized Chinese desert sites and linearly regressed the AVHRR reflectance with cosine SZA to determine the AVHRR degradation. Heidinger et al.²⁰ found a second order relationship between AVHRR TOA reflectance and SZA over Libya-4. Wu et al.³¹ derived a second order cosine SZA relationship with TOA reflectance over Libya-4 and found that the observed Terra and Aqua-MODIS $0.64\text{-}\mu \mathrm{m}$ BRDF corrected reflectance record had a trend standard error of 1.0% and 0.9%, respectively. Initially, the CERES IGCG simply regressed near-nadir observed TOA radiances by cosine SZA (CSZA), but separately for forward and backscatter conditions, because the azimuthally dependent radiances diverged especially for shorter wavelengths.^32–34 More recently, the CERES IGCG began binning the near-nadir observations by repeat day or by scan angle.^35,36

Although the Libya-4 atmospheric column is extremely dry, water vapor absorption impacts many of the near infrared (NIR) and short-wave infrared (SWIR) band TOA radiances, especially during the summer when hot temperatures allow for greater water vapor burden even for very low relative humidity conditions.^20,28 The Libya-1 and Libya-4 PICS have the least seasonal precipitable water (PW) and ozone (O3) variation of all Saharan desert PICS (Bacour et al.¹⁹ see their Fig. 10). Yu and Wu³⁷ noted that the $11\text{-}\mu \mathrm{m}$ minus $12\text{-}\mu \mathrm{m}$ brightness temperature difference was correlated with PW as well as the AVHRR $0.86\text{-}\mu \mathrm{m}$ band TOA reflectance. The PW-corrected AVHRR $0.86\text{-}\mu \mathrm{m}$ band reflectance uncertainty was reduced from $\sim 3\%$ to 1.7%. Multiple studies utilize Hyperion spectra to account for atmospheric absorption and spectral response function differences dependent on viewing and solar angles.^29,38–41 Other studies employed atmospheric correction to account for the ozone and water vapor absorption as well as aerosol scattering over Libya-4^42–45 and Dunhuang⁴⁶ deserts. The CERES IGCG has recently correlated the RSB water vapor absorption with multiple PW sources.^35,36

1.2.

Outline

This work is a continuation of the previous Libya-4 characterization studies of Doelling et al.³⁵ performed to facilitate the CERES IGCG objectives of MODIS and VIIRS imager stability assessment and radiometric scaling. The Terra/Aqua MODIS and NPP/NOAA20 VIIRS reflective solar bands (RSB) are characterized across all angular conditions observed during the 16-day repeat cycle. The previous work characterized only near nadir measurements having view angles less than 30 deg. The utilization of the full scan offers daily observations over the site, thereby greatly increasing the sampling and allowing for a more rigorous identification of pristine clear-sky events. The clear-sky identification is based on the application of both broad and dynamic angular dependent spatial homogeneity thresholds, applied to both the $0.65\text{-}\mu \mathrm{m}$ and $1.61\text{-}\mu \mathrm{m}$ channels, which are the channels least impacted by the atmospheric column. The previous study simply utilized a $0.65\text{-}\mu \mathrm{m}$ channel broad spatial homogeneity threshold to determine clear-sky conditions. This study applies the same characterization across the MODIS and VIIRS imagers, whereas the previous study simply focused on Aqua-MODIS.

Each angular bin, defined by viewing angle and by forward or backscatter azimuth conditions, is characterized by a second order cosine solar zenith relationship, with inclusion of linear PW, O3, and aerosol optical depth (AOD) atmosphere correction terms. This study greatly reduced the number of terms required for characterization compared with the previous effort and thereby mitigates the risk of overfitting. This study also linearly correlated several PW, O3, and AOD daily measurement datasets with the observed TOA channel radiance to find the datasets providing the most effective atmospheric corrections. The atmospheric correction is further validated by the consistency of the relative reflectance trending across RSB and scan angle. Figure 1 provides a flowchart of both the overall methodology and clear-sky identification. The MODIS and VIIRS analogous channel Libya-4 normalized reflectance temporal trends are also compared. Finally, a table of MODIS and VIIRS RSB trend standard errors, with and without atmospheric correction, are provided for full scan and near nadir observations. Section 2 describes the Libya-4 methodology, and Sec. 3 gives the validation results. Concluding remarks and future research are discussed in Sec. 4.

Fig. 1

Schematic of the methodology described in Sec. 2 and the clear-sky identification denoted in Sec. 2.3.

2.1.

Data

The Terra-MODIS, Aqua-MODIS, NPP-VIIRS, and NOAA20-VIIRS data used in this study are a reformatted subset of the original NASA MODIS C6.1 L1b product and the NASA VIIRS Land Science Investigator-led Processing System (Land SIPS) L1b Collection 2. In support of CERES, these data are generated at the Atmospheric Science Data Center Distributed Active Archive Center (ASDC DAAC) located at the NASA Langley Research Center. The MODIS C6.1 radiances were converted to C7 using LUTs provided by the MODIS characterization support team (MCST) and are a function of scan angle, mirror side, and day. The differences between C7 and C6.1 calibrations are outlined in previous studies.⁴⁷ The CERES project subsamples the MODIS 1-km and VIIRS L1B 750-m pixels by utilizing every alternate line and pixel. Table 1 describes the MODIS and VIIRS channel naming convention referred to in the text.

Table 1

MODIS and VIIRS RSB nomenclature as a function of wavelength.

This study compares several sources of PW, O3, and AOD measurements utilized in the Libya-4 PICS atmospheric correction strategy. The CERES project utilizes the Goddard Space Flight Center (GSFC) Global Modeling and Assimilation Office (GMAO) Goddard Earth Observing System 5.4.1 reanalysis product to describe the regional hourly atmospheric profile.⁴⁸ The GMAO PW (PW-GMAO) is based on the assimilation of Aqua-atmospheric infrared sounder and the advanced microwave sounding unit-A observations, as well as ground-based soundings. The GMAO ozone (O3-GMAO) is derived from NPP and NOAA20 - ozone mapping and profiler suite and Aura-ozone monitoring instrument UV observations. The GMAO assimilations strive for a consistent 40-year reanalysis dataset incorporating overlapping satellite records.

The daily Aqua-MODIS MYD08-D3 product provides PW, ozone, and AOD retrievals. The MYD08-D3 NIR approach (PW-RSB) compares the 0.905, 0.936, and $0.94\text{-}\mu \mathrm{m}$ water vapor absorbing channels with the 0.865 and $1.24\text{-}\mu \mathrm{m}$ atmospheric window channels.^49,50 The MYD08-D3 IR (PW-IR) approach mainly relates the 6.72 and $7.33\text{-}\mu \mathrm{m}$ water vapor channels with the $8.55\text{-}\mu \mathrm{m}$ window channel.⁵¹ The MYD08-D3 ozone (O3-MYD08) is retrieved from the MODIS $9.6\text{-}\mu \mathrm{m}$ channel.⁵¹ The NASA-GSFC Aqua-MODIS MYD04_L2 product combines the dark target⁵² and deep blue⁵³ AOD obtained during clear-sky events (AOD-MYD08). Over Libya-4, the MODIS deep blue algorithm is employed. The multi-scale atmospheric transport and chemistry (MATCH) Eulerian three-dimensional model assimilates MODIS AOD retrievals and advects the aerosols in 3-hourly time-steps (AOD-MATCH).⁵⁴ The desert dust-sizes are determined by wind speed. The AOD-MATCH, PW-GMAO, and O3-GMAO measurements are obtained from the CERES SYN1deg product⁵⁵ whereas the PW-RSB, PW-IR, O3-MYD08, and AOD-MYD08 are obtained from the MYD08-D3 product⁵⁶ measurements.

A comparison of the PW measurements over Libya-4 is shown in Fig. 2. The daily PW-GMAO and PW-RSB clear-sky identified (see Sec. 2.3) observations are linearly correlated, revealing that the PW-RSB values are increasing at a greater rate than the PW-IR values [Fig. 2(a)]. PW-GMAO seems to have a nonlinear relationship and overall lower values than the PW-RSB values [Fig. 2(b)], but higher values than those of PW-IR [Fig. 2(c)]. The PW dataset that provides the greatest correlation with the daily observed clear-sky identified TOA radiances is discussed in Sec. 2.5.

Fig. 2

Comparison of the Libya-4 daily (a) PW-RSB and PW-IR, (b) PW-RSB and PW-GMAO, and (c) PW-IR and PW-GMAO clear-sky identified observations over the NPP record.

2.2.

Libya-4 Angular Bins

The CERES IGCG Libya-4 PICS methodology follows Bhatt et al.³³ and Doelling et al.³⁵ The instantaneous Libya-4 MODIS and VIIRS pixel radiances are averaged over the $1\mathrm{deg}\times 1\mathrm{deg}$ region of interest (ROI) centered at 28.6°N latitude and 23.4°E longitude. If the L1b granule boundaries intersect the ROI, the two granules are first concatenated to ensure complete spatial sampling across the entire ROI. The instantaneous observations are easily stratified into angular bins by view zenith angle (VZA) and relative azimuth angle (RAA), as shown in Fig. 3. See also Fig. 1 in Chang et al.²⁸ The Terra and Aqua orbit altitude of 720-km limits the MODIS Earth view scan angle to 65 deg [Fig. 3(a)]. There are eight backscatter VZA and eight forward scatter VZA bins, with each bin representing one repeat cycle day. Due to the higher 833-km altitude of the NPP and NOAA20 satellites, the VIIRS swath encompasses a scan angle of 70 deg and provides 21 angular bins, 16 for each repeat cycle day and 5 from neighboring orbits where the swaths overlap Libya-4.

Fig. 3

Libya-4 (a) Aqua-MODIS and (b) NPP-VIIRS VZA and RAA, with RAA of 0 deg representing forward scatter, observed frequency distributions. The MODIS and VIIRS observations are stratified into 16 and 21 angular bins, respectively.

2.3.

Libya-4 Clear-Sky Identification

The key to successful clear-sky identification over Libya-4 is finding RSBs that are the least susceptible to gaseous absorption and aerosol scattering over the Libya-4 target, especially when observing oblique angular bins. The 0.86, 1.24, and $2.15\text{-}\mu \mathrm{m}$ channels are impacted by water vapor absorption, whereas the 0.48 and $0.55\text{-}\mu \mathrm{m}$ channels are impacted by aerosols and ozone absorption, respectively. This leaves the 0.65 and $1.6\text{-}\mu \mathrm{m}$ channels to determine clear-sky conditions. This study relies on spatial homogeneity filters as an effective way to identify clear-sky events.⁵⁷ Clear-sky events are more spatially uniform than cloudy events. Also, the spatial clear-sky uniformity is assumed to remain consistent over time. The clear-sky events in Fig. 4 are revealed by the densely clustered daily homogeneity values, whereas the cloudy events have a wide range of homogeneity values, which exceed the clear-sky values. Two spatial homogeneity filters are evaluated: the pixel radiance standard deviation over the ROI (${\text{spatial}}_{\mathrm{SD}}$) and a spatial homogeneity factor computed by dividing the ${\text{spatial}}_{\mathrm{SD}}$ by the mean radiance (${\text{spatial}}_{\mathrm{HOM}}$). It was found for channels with wavelengths less than $1\mu \mathrm{m}$ that the ${\text{spatial}}_{\mathrm{SD}}$ has a distinct seasonal fluctuation, with greater values in winter and smaller values in summer. A simple ${\text{spatial}}_{\mathrm{SD}}$ threshold determined for winter conditions would not filter out all cloudy events during the summer. For visible wavelengths ($<1\mu \mathrm{m}$), the ${\text{spatial}}_{\mathrm{HOM}}$ was applied. However, the ${\text{spatial}}_{\mathrm{SD}}$ for wavelengths greater than $1\mu \mathrm{m}$ had a reduced seasonal oscillation, and for these wavelengths, a ${\text{spatial}}_{\mathrm{SD}}$ filter was applied.

Fig. 4

Libya-4 NPP-VIIRS daily observed (a) M3 ${\text{spatial}}_{\mathrm{HOM}}$, (b) M5 ${\text{spatial}}_{\mathrm{HOM}}$, (c) M10 ${\text{spatial}}_{\mathrm{SD}}$, and (d) M11 ${\text{spatial}}_{\mathrm{SD}}$. Panels (e)–(h) same as panels (a)–(d), except after application of the broad and dynamic clear-sky filters. The clear-sky spatial homogeneity values are much lower and similar (densely packed) in time than the cloudy events. The broad homogeneity filter threshold is overlaid in panels (b) and (c).

The Libya-4 clear-sky filtering relies on a two-step approach that consists of broad and dynamic clear-sky filtering. The following broad clear-sky thresholds are utilized; $0.65\text{-}\mu \mathrm{m}$ (VIIRS M5) ${\text{spatial}}_{\mathrm{HOM}}<0.035$, $1.61\text{-}\mu \mathrm{m}$ (VIIRS M10) ${\text{spatial}}_{\mathrm{SD}}<1.0{\mathrm{Wm}}^{-2}{\mathrm{sr}}^{-1}{\mu \mathrm{m}}^{-1}$, $11.0\text{-}\mu \mathrm{m}$ (VIIRS M15) ${\text{spatial}}_{\mathrm{SD}}<2\mathrm{K}$, AOD-MYD08 $<0.6$, AOD-GMAO $<0.75$, and O3-MYDO $<300$ DU. The IR ${\text{spatial}}_{\mathrm{SD}}$ was more effective than the visible ${\text{spatial}}_{\mathrm{SD}}$ at detecting high thin cirrus and their associated shadows. The AOD and O3 parameter thresholds were necessary for removing retrieval outliers, which otherwise may impact the subsequent atmospheric corrections. Figures 4(b) and 4(c) reveal the NPP-VIIRS daily M5 ${\text{spatial}}_{\mathrm{HOM}}$ and M10 ${\text{spatial}}_{\mathrm{SD}}$ observations overlaid by their associated homogeneity thresholds. Note that the NPP-VIIRS M3 ($0.49\text{-}\mu \mathrm{m}$) channel has a much larger clear-sky ${\text{spatial}}_{\mathrm{HOM}}$ range [Fig. 4(a)] than that of the M5 channel [Fig. 4(b)] because it is impacted by both aerosols and Rayleigh scattering. The M3 channel has increased Rayleigh scattering and darker surface reflectance than the M5 channel. Similarly, the M11 channel has a broader clear-sky standard deviation range [Fig. 4(d)] than the M10 channel [Fig. 4(c)] because it is impacted by WV absorption. Overall, the broad clear-sky filters remove $\sim 40\%$ of the observations.

To remove any remaining cloudy events, a dynamic clear-sky filter is applied individually to each angular bin. The angular bin M5 ${\text{spatial}}_{\mathrm{HOM}}$ or M10 ${\text{spatial}}_{\mathrm{SD}}$ mean and standard deviation are computed, after which any observations that exceed twice the standard deviation from the mean are removed. The dynamic filtering process resulted in the removal of an additional 9% of the total observations, predominantly eliminating the winter M5 ${\text{spatial}}_{\mathrm{HOM}}$ spikes [compare Fig. 4(b) with Fig. 4(f) during 2020]. The 0.65 and $1.61\text{-}\mu \mathrm{m}$ dynamic filtering effectively removed most of the daily ${\text{spatial}}_{\mathrm{HOM}}$ outliers across all RSB channels (see Fig. 4, second row). The few remaining M3, M4, M5, and M10 outliers (1.7% of total) are removed utilizing a normalized radiance outlier ($3\sigma$) filter, leaving 49% of the daily observations for stability monitoring.

2.4.

Libya-4 Angular Radiance Models

The MODIS and VIIRS L1b radiances were converted to 1-AU by accounting for the Earth-Sun distance. For each MODIS and VIIRS angular bin, the Libya-4 observed radiances are characterized by applying a second order fit with respect to CSZA utilizing Eq. (1) to obtain the TOA BRDF coefficients ($g_0$, $g_1$, and $g_2$) and thereby compute the predicted radiances (${\mathrm{Rad}}_{\mathrm{pred}}$):

To demonstrate the efficacy of the spatial homogeneity filter for identifying pristine clear-sky conditions, all clear and cloudy Libya-4 daily radiance observations for the most oblique backscatter bin are plotted in Fig. 5 as a function of CSZA and are color coded by their ${\text{spatial}}_{\mathrm{HOM}}$ values. For the M3 and M5 bands shown in Figs. 5(a) and 5(b), respectively, the brighter radiances are associated with greater ${\text{spatial}}_{\mathrm{HOM}}$ (dark red points) than the darker clear-sky radiances, which are highly correlated with CSZA having smaller ${\text{spatial}}_{\mathrm{HOM}}$ (blue points). For the M10 and M11 bands, the clouds are darker than the bright underlying desert. Figures 5(c) and 5(d) show that cloudy M10 and M11 radiances are darker and have larger ${\text{spatial}}_{\mathrm{HOM}}$ (dark red points) than the brighter clear-sky radiances, which have a smaller ${\text{spatial}}_{\mathrm{HOM}}$ (blue points). The second order regression standard error [Eq. (1)] for all observations for this angular bin is shown in Table 2 under the no filter row. The broad filter thresholds (Sec. 2.3) remove most of the cloudy events. The M5 and M10 channel clear-sky dynamic filter reduces the fit standard error below 1%. If the large ${\text{spatial}}_{\mathrm{HOM}}$ events were not related with cloudy events, then the ${\text{spatial}}_{\mathrm{HOM}}$ filtering would not reduce the radiance fit standard error. The success of spatial homogeneity filters avoids utilizing RSB reflectance thresholds, which would otherwise be influenced by wavelength-dependent cloud brightness.

Fig. 5

Libya-4 NPP-VIIRS daily all-sky (not clear-sky filtered) observed (a) M3, (b) M5, (c) M10, and (d) M11 radiances (${\mathrm{Wm}}^{-2}{\mathrm{sr}}^{-1}{\mu \mathrm{m}}^{-1}$) as a function of cos(SZA) for the $68\mathrm{deg}<\mathrm{VZA}<70\mathrm{deg}$ and $90\mathrm{deg}<\mathrm{RAA}<180\mathrm{deg}$ backscatter angular bin. The daily radiances are color coded by the ${\text{spatial}}_{\mathrm{HOM}}$ factor. Note that, for all bands, the cloudy erratic radiances (dark brown points) are easily identified by their large ${\text{spatial}}_{\mathrm{HOM}}$, leaving the clear-sky radiances, which are well correlated.

Table 2

Libya-4 NPP-VIIRS Eq. (1) regression standard error in % for the 68  deg<VZA<70  deg and 90  deg<RAA<180  deg backscatter angular bin shown in Fig. 3. The no filter regression corresponds to the Fig. 5 plots, which utilizes all daily (all-sky) observations. The broad filter applies the static spatial homogeneity thresholds, whereas clear sky includes both the static and dynamic spatial homogeneity thresholds. Note that, for non-absorbing bands (M5 and M10), the clear-sky fit standard error is less than 1%.

The Libya-4 NPP-VIIRS M3, M5, and M11 individual angular bin [Eq. (1)] models are shown in Fig. 6. The M3 angular bin models reveal that the Libya-4 radiance is dependent on the azimuth and view angle. The M5 shows very similar angular bin models and is the most Lambertian channel. M11 shows similar radiances for small CSZA and diverges significantly for near overhead Sun conditions. The RSB angular bin models that are not shown gradually transition by wavelength between Fig. 6 angular models. It may be possible to group the near nadir angular bin radiances into a singular angular bin model [Eq. (1)]; however, by including large view angle bins, the radiances cannot be resolved with a singular angular bin model. The Roujean semi-empirical model⁵⁸ was compared with the second order CSZA model. For near nadir and large VZA bins, the trend standard errors were comparable; however, for the midrange VZA bins, the Roujean trend standard errors were twice as large as the CSZA model (not shown). The second order CSZA relationship describes the observed clear-sky radiances with the fewest terms.

Fig. 6

Libya-4 NPP-VIIRS daily predicted (a) M3, (b) M5, (c) M11, radiance (${\mathrm{Wm}}^{-2}{\mathrm{sr}}^{-1}{\mu \mathrm{m}}^{-1}$) as a function of CSZA for the 21 angular bins color coded by scan angle; the blue bins are in the backscatter direction or the left side of the scan, the red bins are in the forward scatter direction or right side of the scan, and nadir scans are shown in green.

2.5.

Libya-4 Atmospheric Correction

In this section, we correlate the PW, O3, and AOD retrievals with the TOA observed clear-sky radiances as a function of wavelength. The Libya-4 TOA radiances are modeled by applying a second order fit to the observed radiances with respect to CSZA, similar to Eq. (1), and incorporating the linear coefficient terms for PW, O3, and AOD, as illustrated in Eq. (2). Equation (2) coefficients are computed for each of the 16 MODIS and 21 VIIRS angular bins as

The normalized radiances (${\mathrm{rad}}_{\text{norm}}$) are computed by dividing daily observed clear-sky radiances (${\mathrm{rad}}_{\mathrm{obs}}$) by the predicted or modeled radiance (${\mathrm{rad}}_{\mathrm{mod}}$) from Eq. (2), as

The linear regression standard error (${\text{trend}}_{\mathrm{SE}}$) is computed from the clear-sky daily ${\mathrm{rad}}_{\text{norm}}$ observations across all angular bins over the record. To determine the atmospheric correction impact, Eq. (2) coefficients are computed for a single atmospheric parameter at a time, and the resulting ${\text{trend}}_{\mathrm{SE}}$ is compared with the ${\text{trend}}_{\mathrm{SE}}$ without utilizing atmospheric terms [Eq. (1)]. A comparison of the 3 PW sources of the ${\text{trend}}_{\mathrm{SE}}$ in Table 3 shows that the PW-RSB has the lowest ${\text{trend}}_{\mathrm{SE}}$ for the M7 ($0.86\mu \mathrm{m}$), M8 ($1.24\mu \mathrm{m}$), and M11 ($2.25\mu \mathrm{m}$) water vapor absorption channels. This is not surprising as the RSB channels were used to retrieve the water vapor in the PW-SRB dataset. Note that the M7 and M8 PW-RSB ${\text{trend}}_{\mathrm{SE}}$ was very close to the M10 ${\text{trend}}_{\mathrm{SE}}$, which was least impacted by the atmospheric column and indicates that the PW-RSB dataset was nearly optimal in estimating the radiance impact of PW.

Table 3

Libya-4 NPP-VIIRS Eq. (2) regression standard error in % across all angular bins as a function of a single atmospheric dataset (see Sec. 2.1). The CSZA-only represents the Eq. (1) linear regression standard error (trendSE) in % of the clear-sky daily radnorm observations. The bold text denotes the trendSE reduction in % by adding a single atmospheric correction term in Eq. (2) compared with the CSZA-only trendSE. For comparison, M5 and M10 bands, which are the least sensitive to the atmospheric column, have a CSZA-only trendSE of 0.686% and 0.636%, respectively.

The O3-GMAO UV based ozone retrieval is slightly more effective in predicting the M4 ozone absorption impact than the O3-MYD08 $9.6\text{-}\mu \mathrm{m}$ ozone retrieval. The ozone impact was minimal due to the small seasonal ozone fluctuations over Libya-4. The AOD-MYD08 is slightly more effective in determining the AOD impact on the M3 ($0.49\mu \mathrm{m}$) radiance because the AOD-MYD08 aerosols were directly retrieved during clear-sky events. Note that the AOD and O3 atmospheric corrections were minimal compared with the PW atmospheric correction reduction. Finding a more robust ozone and aerosol retrieval product or utilizing nonlinear terms with TOA reflectance is left for future studies. Similar results were obtained for NOAA20-VIIRS, Terra-MODIS, and Aqua-MODIS. There is no advantage in using the Terra-based PW-RSB, PW-IR, or AOD over their Aqua counterparts.

Based on the analysis in Table 3, the PW-RSB, O3-GMAO, and AOD-MYD08 datasets were used to determine the PW, O3, and AOD terms in Eq. (2). Figure 7(a) shows the Libya-4 NPP VIIRS M11 clear-sky daily ${\mathrm{rad}}_{\text{norm}}$ without applying the atmospheric correction terms, which reveals a large seasonal oscillation of $\pm 5\%$, with the positive ${\mathrm{rad}}_{\text{norm}}$ values occurring in spring and the negative ${\mathrm{rad}}_{\text{norm}}$ values occurring in fall. A slight downward trend is observed (0.1%/year). Figure 7(b) shows the clear-sky daily atmosphere-corrected ${\mathrm{rad}}_{\text{norm}}$ observations. No discernable seasonal oscillation is perceived, and the downward trend is reduced by half. Figure 7(b) reveals that all daily angular bin ${\mathrm{rad}}_{\text{norm}}$ are equally dispersed, indicating that the atmospheric corrections are effective across all angular bins. The atmospheric correction reduced the ${\text{trend}}_{\mathrm{SE}}$ by $\sim 63\%$, which is almost entirely due to the PW contribution for this band. The downward trend is discussed in the following sections and is more than likely due to the Libya-4 surface natural variability.

Fig. 7

Libya-4 NPP-VIIRS M11 clear-sky daily ${\mathrm{rad}}_{\text{norm}}$ observations (a) CSZA-only terms [Eq. (1)] and (b) with atmospheric correction terms [Eq. (2)]. The daily ${\mathrm{rad}}_{\text{norm}}$ observations are color coded by scan angle. The CSZA-only and with atmospheric correction ${\text{trend}}_{\mathrm{SE}}$ are 1.92% and 0.71%, respectively, over the 11-year record. Similarly, the linear fit slopes are $-0.011$ and $-0.0059$, respectively.

3.1.

Libya-4 Scan Angle Consistency

The atmospheric correction is compared across angular bins. The atmospheric column radiance contribution is a function of the gaseous absorption contained within the channel spectral response and the combined SZA and VZA pathlength. If the atmospheric correction were perfect, the ${\mathrm{rad}}_{\text{norm}}$ temporal variability across the observed pathlengths or angular bins would be similar and would follow the Libya-4 surface reflectance natural variability given that surface reflectance is near Lambertian. The seasonal or interannual gaseous absorption or aerosol variability, after atmospheric correction, should not be manifested in the temporal record. To facilitate the angular bin consistency comparison of the observations in Fig. 7, a 365-day running mean is applied to smooth the individual angular bin daily ${\mathrm{rad}}_{\text{norm}}$ observations and is shown in Fig. 8. Figures 7(a) and 7(b) daily ${\mathrm{rad}}_{\text{norm}}$ observations plots can be compared with Figs. 8(c) and 8(f) plots, where the latter plots have the 365-day running mean applied. The atmosphere correction terms reduced the overall temporal variability. Also, Fig. 8 atmospheric corrected plots show in-synch angular bin temporal trending, especially for the 1.24 and $2.25\text{-}\mu \mathrm{m}$ channels. The $0.55\text{-}\mu \mathrm{m}$ ozone correction has tightened some of the noisier angular bin trends [compare Fig. 8(a) with Fig. 8(d)]. Also, the atmospheric-corrected $0.55\text{-}\mu \mathrm{m}$ channel has similar temporal trending as the water vapor channels in Fig. 8.

Fig. 8

Libya-4 NPP-VIIRS (a) M4, (b) M8, (c) M11, 365-day running mean ${\mathrm{rad}}_{\text{norm}}$ observations using CSZA-only terms [Eq. (1)], color coded by angular bin. Panels (d)–(f) same as panels (a)–(c), except including atmospheric correction terms [Eq. (2)].

If the atmospheric radiance contributions were perfectly removed, then the ${\text{trend}}_{\mathrm{SE}}$ value should be similar across all RSBs for a given angular bin. Figure 9(a) shows the NPP-VIIRS RSB ${\text{trend}}_{\mathrm{SE}}$ as a function of angular bin without atmospheric correction. The ${\text{trend}}_{\mathrm{SE}}$ increases with increasing VZA, especially for backscatter conditions. As expected, the M5 and M10 have the lowest ${\text{trend}}_{\mathrm{SE}}$ because these bands are the least impacted by atmosphere, whereas the bands impacted by water vapor and AOD have the highest ${\text{trend}}_{\mathrm{SE}}$. Figure 9(b) shows the NPP-VIIRS band ${\text{trend}}_{\mathrm{SE}}$ with atmospheric correction. The atmospheric correction reduced the M7, M8, and M11 ${\text{trend}}_{\mathrm{SE}}$ and are more in line with the non-absorbing M5 and M10 channels. However, for the M3 and M4 channels, the ${\text{trend}}_{\mathrm{SE}}$ was not significantly reduced with AOD and O3 atmospheric correction terms, suggesting that the AOD and O3 sources exhibited only minimal linear correlation with TOA radiance. The overall Libya-4 ${\text{trend}}_{\mathrm{SE}}$ can be reduced by restricting the number of large VZA bins.

Fig. 9

Libya-4 NPP-VIIRS 365-day averaged spectral band ${\text{trend}}_{\mathrm{SE}}$ using (a) CSZA-only terms [Eq. (1)] and (b) including atmospheric correction terms [Eq. (2)], color coded by spectral band as denoted by the legend. Note that many of the water vapor absorbing band angular bin ${\text{trend}}_{\mathrm{SE}}$ are more similar to the M5 and M10 band ${\text{trend}}_{\mathrm{SE}}$ after applying atmospheric correction terms.

3.2.

Libya-4 RSB Consistency

We now compare the atmospheric correction across the RSB. Without atmospheric corrections, the disparate channel gaseous absorption ${\mathrm{rad}}_{\text{norm}}$ temporal fluctuations would be uncorrelated, if we assume that the Libya-4 surface reflectance temporal variability is minimally dependent on wavelength. Figure 10 left and right column plots display the MODIS and VIIRS RSB ${\mathrm{rad}}_{\text{norm}}$ observations without and with atmospheric corrections, respectively. Note the large VIIRS M11 and MODIS B7 (yellow line) and B17 (gray line) water vapor absorption band ${\mathrm{rad}}_{\text{norm}}$ temporal fluctuations. After atmospheric correction, the large fluctuations are greatly reduced and fluctuate in tandem with the other channels. Except for Terra-MODIS, the RSB ${\mathrm{rad}}_{\text{norm}}$ fluctuations are very consistent, thereby validating that the RSB and angular-bin-specific atmospheric correction coefficients. The Terra-MODIS instrument onboard calibration system anomalies caused the spurious L1b radiance drifts. Furthermore, the Terra-MODIS solar diffuser door has remained in the open position since 2003, which increased the solar diffuser degradation for shorter wavelengths by 20% when compared with Aqua-MODIS. Finally, the 2016 Terra spacecraft anomaly increased the electronic cross-talk, especially for the SWIR bands.⁵⁹

Fig. 10

Libya-4 Terra-MODIS 365-day averaged ${\mathrm{rad}}_{\text{norm}}$ observations using (a) CSZA-only terms [Eq. (1)] and (b), including atmospheric correction terms [Eq. (2)], color coded by the spectral band according to the legend located in the right column plots. Panels (c) and (d) same as panels (a) and (b) except for Aqua-MODIS. Panels (e) and (f) same as panels (a) and (b) except for NPP-VIIRS. Panels (g) and (h) same as panels (a) and (b) except for NOAA20-VIIRS. In comparing the CSZA-only (left column) and with atmospheric correction (right column) plots, note the improved RSB consistency attained by correcting for the atmosphere.

Note that the same positive $+0.5\%$ ${\mathrm{rad}}_{\text{norm}}$ inflection during 2014 is present in the Terra-MODIS [Fig. 10(b)], Aqua-MODIS [Fig. 10(d)], and NPP-VIIRS [Fig. 10(f)] records. It seems that the RSB ${\mathrm{rad}}_{\text{norm}}$ fluctuations are correlated across sensors. The 2014 ${\mathrm{rad}}_{\text{norm}}$ spike also corresponds to the increased solar radiation exposure that degraded the MODIS and VIIRS solar diffusers.⁶⁰ However, the solar radiation exposure event should be accounted for in the L1B product by the onboard calibration system solar diffuser stability monitor.

Here, we verify if the atmospheric-corrected ${\mathrm{rad}}_{\text{norm}}$ fluctuations are comparable across the MODIS and VIIRS imagers. The 0.65 and $0.86\text{-}\mu \mathrm{m}$ MODIS and VIIRS band ${\mathrm{rad}}_{\text{norm}}$ observations are presented in Figs. 11(a) and 11(b), respectively, which show that, indeed, the ${\mathrm{rad}}_{\text{norm}}$ temporal fluctuations are in harmony. However, the $1.64\text{-}\mu \mathrm{m}$ band ${\mathrm{rad}}_{\text{norm}}$ temporal fluctuations [Fig. 11(c)] differ slightly from their visible counterparts before 2012. The 2014 ${\mathrm{rad}}_{\text{norm}}$ inflection is present across all RSB and all sensors, except NOAA20-VIIRS because it was launched in late 2017, suggesting that the atmospheric correction and BRDF inadequacies do not explain the 2014 inflection.

Fig. 11

Libya-4 365-day averaged atmospheric corrected ${\mathrm{rad}}_{\text{norm}}$ observations for the (a) $0.65\text{-}\mu \mathrm{m}$ wavelength color coded by satellite and RSB as denoted in the legend. Panel (b) same as panel (a) except for the $0.86\text{-}\mu \mathrm{m}$ bands. Panel (c) same as panel (a) except for the $1.64\text{-}\mu \mathrm{m}$ bands. The MODIS and VIIRS bands are listed in Table 1.

The fact that ${\mathrm{rad}}_{\text{norm}}$ inflections are in-synch among the MODIS and VIIRS RSB indicates that the inflections accurately describe the Libya-4 site natural variability and that the sensor calibration anomalies are much smaller than the true Libya-4 natural variability. The inflections being similar between RSBs suggests that the Libya-4 natural variability is mainly owed to prevailing winds changing the orientation of sand dunes that define the areal extent of the shadowing.²⁴ This finding suggests that (1) the Libya-4 sand color is not changing in time from dust advection, (2) wet sand evaporates quickly from rain events, and (3) the sand dunes prevent any vegetation from growing.

3.3.

Libya-4 Temporal Stability

We now quantify the Libya-4 stability by the imager, RSB, and atmospheric correction. Table 4 quantifies the MODIS and VIIRS channel-specific ${\text{trend}}_{\mathrm{SE}}$ without (first row of Table 2) and with (second row of Table 2) atmospheric corrections, with the percent ${\text{trend}}_{\mathrm{SE}}$ reduction due to atmospheric correction being given in the third row in red text. The NPP and NOAA20 VIIRS M8 and M11 had a $\sim 40\%$ and $\sim 62\%$ ${\text{trend}}_{\mathrm{SE}}$ reduction by including the PW term. The Terra and Aqua MODIS B2, B17, B5, and B7 ${\text{trend}}_{\mathrm{SE}}$ was reduced by $\sim 35\%$, 80%, 30%, and 67% with PW correction, respectively. The $1.61\text{-}\mu \mathrm{m}$ wavelength, which is least impacted by the atmosphere, has the smallest RSB ${\text{trend}}_{\mathrm{SE}}$, which are within 0.6% and probably represents the Libya-4 natural variability. Most of the water vapor absorbing channel ($>0.8\mu \mathrm{m}$) ${\text{trend}}_{\mathrm{SE}}$ after atmospheric correction are within 0.7%, indicating that PW mostly accounts for the daily TOA radiance variability. The ${\text{trend}}_{\mathrm{SE}}$ atmospheric correction reduction was less than 20% for wavelengths less than $0.6\mu \mathrm{m}$, which implies that a more effective AOD and ozone source or a better pristine clear-sky filtering approach is needed, which is left for future studies. The larger Terra-MODIS $0.48\text{-}\mu \mathrm{m}$ channel ${\text{trend}}_{\mathrm{SE}}$, when compared with Aqua-MODIS, are due to unresolved polarization issues.⁵⁹ The similarity between MODIS and VIIRS ${\text{trend}}_{\mathrm{SE}}$ suggests that the MODIS L1b C6.1 to C7 temporal radiometric adjustments are effective.⁶¹ In general, the Libya-4 ${\text{trend}}_{\mathrm{SE}}$ is between 0.6% and 0.75% for wavelengths greater than $0.6\mu \mathrm{m}$ and between 0.75% and 1.0% for wavelengths less than $0.6\mu \mathrm{m}$.

Table 4

Libya-4 Terra-MODIS, Aqua-MODIS, NPP-VIIRS, and NOAA20-VIIRS trendSE by RSB. For each sensor, the first row is based on the CSZA trendSE, the second row is the trendSE with atmospheric corrections, and the third row (bold text) is the atmospheric correction trendSE reduction in %.

We now perform the Libya-4 stability assessment for viewing angles less than 35 deg. It was noted in Sec. 3.1 and Fig. 9(b) that the NPP-VIIRS near-nadir angles have smaller ${\text{trend}}_{\mathrm{SE}}$. Similar scan angle dependent ${\text{trend}}_{\mathrm{SE}}$ was obtained for MODIS and NOAA20-VIIRS (not shown). Near nadir observations avoid oblique path lengths that require larger atmospheric corrections and three-dimensional sand dune effects. Limiting the viewing angle to 35 deg or less reduces the sampling by $\sim 66\%$, which decreases the number of clear-sky observations over the 16-day repeat cycle, but it is countered by the increased temporal stability of the near nadir observations. Also, the long-term MODIS Libya-4 reflectance trending by scan angle is much noisier for oblique viewing angles (see Fig. 13, Xiong et al.⁵⁹). The VIIRS instrument was designed to mitigate the response versus scan angle variations by utilizing a rotating telescope assembly and a half angle mirror, unlike MODIS, which relies on a scan mirror.⁴⁷

Table 5 shows the Libya-4 ${\text{trend}}_{\mathrm{SE}}$ for viewing angles less than 35 deg and can be directly compared with Table 4, which includes all angles. The RSB dependent atmospheric correction ${\text{trend}}_{\mathrm{SE}}$ percentage reductions are similar for near-nadir and full scan conditions. This result suggests that performing atmospheric corrections is effective not just for oblique angle observations but also for nadir observations. A 5% to 23% ${\text{trend}}_{\mathrm{SE}}$ reduction was observed by limiting the VZA to 35 deg, depending on the channel. The VIIRS non-absorbing $1.61\text{-}\mu \mathrm{m}$ channel ${\text{trend}}_{\mathrm{SE}}$ for VZA less than 35 deg is within 0.50%, whereas for all viewing angles, it was within 0.60%, indicating that the surface reflectance variability is slightly larger for oblique views. Overall, the Libya-4 ${\text{trend}}_{\mathrm{SE}}$ for VZA less than 35 deg was between 0.5% and 0.8%, except for the Terra-MODIS B3 channel.

Table 5

Same as Table 4 except for angular bins with VZA<35  deg.