Comparative analysis of GOME and SCIAMACHY reflectance over Pseudo-Invariant Calibration Sites: implications for spectrometers cross-calibration

Accurate radiometric cross-calibration is essential for ensuring the consistency and interoperability of multi-sensor satellite observations. Vicarious calibration, a widely adopted approach, relies on the temporal stability of desert-based Pseudo-Invariant Calibration Sites (PICS). However, these sites are limited in spatial extent and have not been systematically examined for cross-calibration with larger pixels, whose dimensions exceed those of PICS. This study established a statistical framework to advance cross-calibration of spectrometers over PICS and their surrounding areas. The framework included performance comparisons of different satellite instruments and the identification of reference and constrained sensors. The temporal stability of PICS across various spectral bands was reassessed using observations from the reference sensor. Stability scores (SS) were derived from a combination of statistical indicators designed to capture temporal variability, distribution symmetry, the occurrence of anomalies, and long-term shifts in the observations.

Decades of top-of-atmosphere (TOA) reflectance data in the ultraviolet, visible, and near-infrared (UV/VIS/NIR) ranges, collected by the Global Ozone Monitoring Experiment (GOME) and Scanning Imaging Absorption Spectrometers for Atmospheric Chartography (SCIAMACHY) spectrometers over 20 PICS sites, were analyzed. The results revealed significant degradation in GOME, particularly in the UV band, during 2001 and toward the end of its mission, as evidenced by positively skewed, heavy-tailed distributions. In contrast, SCIAMACHY observations were more uniform and stationary, indicating their greater suitability for assessing PICS stability. The investigated sites were ranked based on the average of SS across the entire investigated spectrum as well as for each spectral band. The analysis revealed a wide range of SS among PICS, including intra-site variations across spectral bands. While some sites demonstrated consistently high stability, many were found to be 2–3 times less stable than the most stable sites. Among these robust sites, some have not been recommended by the calibration communities and should be given further consideration. In addition, the stability of the sites based on TOA reflectance showed a clear spectral band dependency, with enhanced stability in the NIR relative to the UV and VIS bands. These findings underscore the importance of regularly evaluating PICS and of considering spectral-band-specific performance when selecting calibration sites and performing the cross-calibration of spectrometers over PICS.

1 Introduction

In recent decades, the number of satellite missions has significantly increased, leading to a paradigm shift across various applications. This expansion has resulted in an abundance of remote sensing data and a rapid growth in data archives. Long-term monitoring applications, such as those focused on atmospheric gases, require continuous time series data that often extend beyond the operational lifespan of individual satellite missions (Cabor et al., 1998). Ensuring radiometric consistency among these datasets is critical for integrating observations from different satellite missions. Furthermore, it is fundamental for developing high data quality to perform quantitative remote sensing (Helder et al., 2013).

Various calibration strategies can be used to ensure consistency across both temporal and spectral domains, ranging from pre-flight calibration (Lichtenberg et al., 2006; Coldewey-Egbers et al., 2008), in-flight calibration (Pagano et al., 2003; Ono et al., 1996), and vicarious calibration (Cabor et al., 1998; Wu et al., 2020; Helder et al., 2013; Chander et al., 2013; Ye et al., 2021). Notably, vicarious calibration is independent of the onboard calibration systems and provides an external means of validating and correcting satellite measurements.

Pseudo-invariant Calibration Sites (PICS) are among the most widely utilized targets in vicarious calibration methods (Kabir et al., 2020; Hasan et al., 2019), and absolute radiometric calibration (Helder et al., 2013). A set of 20 desert sites, each measuring 100×100 km², was proposed across Africa and the Arabian Peninsula by Cosnefroy et al. (1996). Several of these sites have been recommended as reference calibration and validation test sites by the Committee on Earth Observation Satellites (CEOS), specifically by the Working Group on Calibration and Validation (WGCV) and the Infrared and Visible Optical Sensors (IVOS) (Thenkabail, 2015; Liang, 2017). As reported in Cosnefroy et al. (1996), the reflectance across each site was spatially uniform, with relative variations better than 3 %. In contrast, the temporal analysis based on cloud-free multitemporal imagery revealed peak-to-peak reflectance variations of approximately 8 % to 15 %.

PICS have been the subject of several recent studies (Khadka et al., 2021; Tuli et al., 2019; Fajardo Rueda et al., 2021; Hasan et al., 2019). These investigations confirmed that PICS remain relevant and reliable for radiometric calibration of optical remote sensing systems. Therefore, any radiometric deviations in the acquired scenes that exceed the expected temporal stability of PICS may suggest potential instrumental anomalies and necessitate further calibration or investigation – this principle forms the basis for employing PICS in radiometric calibration, sensor performance monitoring, and cross-calibration between different sensors.

While PICS are well-established in the radiometric and cross-calibration of imagery sensors, their applications to spectrometer calibration remain largely unexplored. Spectrometer footprints typically have a much larger spatial extent than that of PICS. Therefore, the temporal stability of areas larger than the PICS domain should be assessed before their use in cross-calibration. To address this gap, this research paper aims to establish a foundation for spectrometer cross-calibration through two key objectives. The first objective is to develop a statistical framework for the comparative analysis of spectrometers, to monitor their performance throughout the mission lifetime, and to identify both the reference and the constrained sensor. The second objective is to assess the temporal stability of PICS and surrounding areas across the full spectrum.

In this study, top-of-atmosphere (TOA) reflectance observations from the Global Ozone Monitoring Experiment (GOME) and the Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY) over selected PICS were analyzed. A decade of measurements was used to construct reflectance time series in three spectral regions: ultraviolet (UV), visible (VIS), and near-infrared (NIR). These time series form the basis for evaluating both sensor performance and site stability. Although the derived comparison metrics inherently include contributions from surface properties and atmospheric effects, the analysis focuses on relative temporal variability rather than absolute reflectance levels.

Relative comparisons are particularly suitable over homogeneous desert regions, which are typically remote from coastlines and major emission sources. Under the assumption that atmospheric conditions are broadly comparable across the selected sites, variations in TOA reflectance are expected to be primarily governed by surface characteristics. Consequently, differences in temporal variability can be attributed mainly to site-specific properties or sensor behavior, enabling the identification of the most stable sites for cross-calibration.

Within this framework, the study is guided by two main hypotheses. First, atmospheric contributions over the selected desert calibration sites are assumed to be temporally and spatially comparable. Under stable sensor operation, variability metrics should therefore be similar across sites. Second, the PICS themselves are expected to exhibit consistent temporal behavior due to their comparable environmental and meteorological settings, regardless of the observing sensor. Systematic deviations from this expected consistency may indicate site-specific effects or sensor-related instabilities.

The paper is structured as follows. Section 2 describes the sensors and their data used in this investigation. Section 3 presents the geographical locations of the investigated PICS. The methods are presented in Sect. 4. Results and discussions are presented in Sects. 5 and 6. Finally, the conclusion and outlook in Sect. 7.

2 Sensor and data descriptions

2.1 GOME

GOME was a passive nadir-scanning UV/VIS spectrometer onboard ERS-2 (ESA, 2025b). It provided observations of ozone, NO₂, clouds, aerosols, and other atmospheric parameters with global coverage every 3 d. GOME had forward and backward scanning sequences with 4 pixels per scan. Forward scans included east, nadir, and west pixels (Burrows et al., 1999).

2.2 SCIAMACHY

SCIAMACHY was a passive UV/VIS/NIR imaging spectrometer on ENVISAT (ESA, 2025c). It featured nadir, limb, and sun/moon occultation viewing. Limb and occultation modes were excluded in this study. SCIAMACHY retrieved trace gases, aerosols, radiation, and cloud properties (Bovensmann et al., 1999).

Further details on the characteristics information of each sensor, including the mission life, orbit, and spatial and spectral resolution, are illustrated in Table 1.

Table 1 Comparison of GOME and SCIAMACHY satellite instruments.

2.3 GOME and SCIAMACHY data

Fully calibrated Level-1 products from GOME (Level-1b version 5.2) (ESA, 2025b) and SCIAMACHY (Level-1c version 10) (ESA, 2025c) were used to derive TOA reflectance time series. Both GOME Level-1b and SCIAMACHY Level-1c datasets contained solar irradiance and radiance measurements required for reflectance computation. While cloud parameters were included in the GOME Level-1b files, they were not available in the SCIAMACHY Level-1c products. Consequently, cloud parameters for SCIAMACHY were obtained from the corresponding Level-2 products.

The following is the main data from both sensors used:

• Reflectance: this parameter refers to the fraction of incoming radiation reflected toward the instruments. This consists of surface reflectance and atmospheric scattering. Reflectance serves as a fundamental parameter for the comparison between sensors and the study of the temporal variability of PICS. In this study, the term “reflectance” refers to TOA reflectance, unless stated otherwise.

• Cloud fraction (CF): this parameter ranges from 0 to 1, where 0 indicates clear skies (no clouds) and 1 indicates full cloud cover. The cloud fraction is used to filter out data with a higher CF in the reflectance time series.

• Viewing geometry of spectrometers: the solar zenith angle (SZA) and viewing zenith angle (VZA) are geometric parameters that may influence the temporal variability of reflectance. SZA is the angle between the local vertical and the sun, and VZA is the angle between the local vertical and the satellite line of sight. The effects of different SZA and VZA were investigated for SCIAMACHY reflectance before computation SS in Sect. 5.4.

3 Study areas description

The PICS proposed by Cosnefroy et al. (1996) are shown in Fig. 1. These sites, located primarily in desert regions, were chosen due to their minimal vegetation cover and stable surface reflectance characteristics. These features made them particularly suitable for satellite sensor intercomparison and cross-calibration. Further details on the geographical locations of the PICS are provided in Table S1 in the Supplement.

Figure 1 PICS in Africa and the Arabian Peninsula with the dimensions of 100×100 km². The sites with green dots were endorsed by the CEOS for cross-calibration purposes (CEOS, 2026).

4 Methods

This section outlines the methodological workflows, addressing the two key objectives of the study: (i) comparative assessment of GOME and SCIAMACHY across the PICS and investigated spectral ranges, and (ii) modeling stability score (SS) for each PICS based on derived metric parameters.

4.1 Comparative assessment of GOME and SCIAMACHY

4.1.1 Computation of the reflectance

The reflectance measurements from both sensors were computed with the assumption of a Lambertian surface. The reflectance from GOME (R_g) and SCIAMACHY (R_s) were computed based on Eq. (1):

$$\begin{array}{}\text{(1)}& R\left(\mathit{\lambda }\right)={\displaystyle \frac{\mathit{\pi }\times I_{\text{earth}}\left(\mathit{\lambda }\right)}{\cos \left({\mathit{\varphi }}_o\right)\times E_{\text{sun}}\left(\mathit{\lambda }\right)}}\end{array}$$

where R(λ) is the reflectance as a function of wavelength, I_earth(λ) is the calibrated radiance as measured by the instrument as a function of wavelength, ϕ_o is the solar zenith angle, and E_sun is the solar irradiance as measured by the instrument.

4.1.2 Retrieving of time series of reflectance from GOME and SCIAMACHY data

The reflectance data from the complete archival record of ESA's GOME and SCIAMACHY spectrometers over each PICS were retrieved. The retrieved reflectance data spanned spectral ranges in UV/VIS/NIR. The selected wavelength ranges and the corresponding number of spectral channels are presented in Table 2. In this paper, the term “spectral channel” referred to an individual spectral pixel of the detector, whereas “spectral band” denoted the detector channel. The selected wavelength ranges were identical to those used in the Fundamental Data Record for Atmospheric Composition (FDR4ATMOS) project for performing cross-calibration between GOME and SCIAMACHY (ESA, 2025a; Owda et al., 2026). These spectral ranges contain absorption features of numerous atmospheric gases.

Table 2 Characteristics of the selected wavelength ranges for each spectral band, including the number of spectral channels corresponding to each band for the GOME and SCIAMACHY sensors.

For each satellite overpass, reflectance values were extracted within a 1.5°×1.5° bounding box centered over the PICS. This area was chosen to be larger than the actual PICS extent because the GOME pixel size exceeded the PICS dimensions, resulting in only a single reflectance value per overpass. In contrast, SCIAMACHY could provide multiple reflectance values within the PICS, and their average was used for analysis.

We applied a threshold of CF≤0.25 to filter the time series of both sensors. This threshold was chosen to minimize the influence of cloud contamination on the observed reflectance values, as high CF could significantly alter the magnitude and variability of the measurements. As a result of applying the CF filtering, the number of samples was reduced by an average of 18 % and 8 % for GOME and SCIAMACHY, respectively. The original and filtered sample counts are reported in Table S1.

4.1.3 Deriving temporal stability metrics

For each investigated PICS site (S_k) and wavelength (λ), we extracted a set of statistical and signal-based features that described the temporal variability of the reflectance time series. The extracted features offered a multidimensional representation of the reflectance dynamic for monitoring sensors' performance and identifying the stability of spectral channels and sites. The parameters were categorized into two different groups as follows:

1. Basic statistical parameters. These parameters include the mean (μ), the standard deviation (σ), the coefficient of variation (CV), and the interquartile range (IQR).

   μ describes the average reflectance over time (Eq. 2), σ is used to describe the absolute temporal variability around the mean (Eq. 3), CV describes the degree of variability in the reflectance time series (Eq. 4). IQR measures the statistical dispersion in the middle range of 50 % of the dataset (Eq. 5).

   For each site S_k and λ from a sensor 𝒮 at time t_n, the reflectance time series was denoted as $R_{\mathcal{S},\mathit{\lambda }}^{S_k}\left(t_n\right)$

   $$\begin{array}{}\text{(2)}& {\displaystyle }{\mathit{\mu }}_{\mathcal{S},\mathit{\lambda }}^{S_k}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{n=\mathrm{1}}^N{R}_{\mathcal{S},\mathit{\lambda }}^{S_k}\left(t_n\right)\text{(3)}& {\displaystyle }{\mathit{\sigma }}_{\mathcal{S},\mathit{\lambda }}^{S_k}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{n=\mathrm{1}}^N{\left(R_{\mathcal{S},\mathit{\lambda }}^{S_k}\left(t_n\right)-{\mathit{\mu }}_{\mathcal{S},\mathit{\lambda }}^{S_k}\right)}^{\mathrm{2}}}\text{(4)}& {\displaystyle }{\text{CV}}_{\mathcal{S},\mathit{\lambda }}^{S_k}={\displaystyle \frac{{\mathit{\sigma }}_{\mathcal{S},\mathit{\lambda }}^{S_k}}{{\mathit{\mu }}_{\mathcal{S},\mathit{\lambda }}^{S_k}}}\text{(5)}& {\displaystyle }{\text{IQR}}_{\mathcal{S},\mathit{\lambda }}^{S_k}=Q_{\mathrm{3}}\left(R_{\mathcal{S},\mathit{\lambda }}^{S_k}\right)-Q_{\mathrm{1}}\left(R_{\mathcal{S},\mathit{\lambda }}^{S_k}\right)\end{array}$$

   where N is the number of observations for each reflectance time series. The number of observations for each site is tabulated in Table S1. Q₃ is the 75th percentile, and Q₁ is the 25th percentile.

2. Trend and shape parameters. These parameters include the slope (m), the skewness (${{\mathit{\gamma }}_{\mathrm{1}}}_{\mathcal{S},\mathit{\lambda }}^{S_k}$), and the kurtosis (${{\mathit{\gamma }}_{\mathrm{2}}}_{\mathcal{S},\mathit{\lambda }}^{S_k}$) of the reflectance time series. The slope was estimated by fitting the reflectance to the linear regression model (Eqs. 6 and 7). It describes the long-term changes and trends in the reflectance over time. γ₁ detects asymmetry in the reflectance distribution and gives insights about the anomaly (Eq. 8). γ₂ measures the tailedness of the distribution (Eq. 9).

   $$\begin{array}{}\text{(6)}& {\displaystyle }{R}_{\mathcal{S},\mathit{\lambda }}^{S_k}\left(t_n\right)=a_{\mathcal{S},\mathit{\lambda }}^{S_k}{t}_n+b_{\mathcal{S},\mathit{\lambda }}^{S_k}+{\mathit{\epsilon }}_n\text{(7)}& {\displaystyle }m=a_{\mathcal{S},\mathit{\lambda }}^{S_k}\text{(8)}& {\displaystyle }{{\mathit{\gamma }}_{\mathrm{1}}}_{\mathcal{S},\mathit{\lambda }}^{S_k}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{n=\mathrm{1}}^N{\left({\displaystyle \frac{R_{\mathcal{S},\mathit{\lambda }}^{S_k}\left(t_n\right)-{\mathit{\mu }}_{\mathcal{S},\mathit{\lambda }}^{S_k}}{{\mathit{\sigma }}_{\mathcal{S},\mathit{\lambda }}^{S_k}}}\right)}^{\mathrm{3}}\text{(9)}& {\displaystyle }{{\mathit{\gamma }}_{\mathrm{2}}}_{\mathcal{S},\mathit{\lambda }}^{S_k}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{n=\mathrm{1}}^N{\left({\displaystyle \frac{R_{\mathcal{S},\mathit{\lambda }}^{S_k}\left(t_n\right)-{\mathit{\mu }}_{\mathcal{S},\mathit{\lambda }}^{S_k}}{{\mathit{\sigma }}_{\mathcal{S},\mathit{\lambda }}^{S_k}}}\right)}^{\mathrm{4}}\end{array}$$

4.2 Computing stability score

4.2.1 Refinement reflectance observations

To ensure SS reliability represents the temporal variability of the sites, the following steps were applied:

• Exclusion of the O₂ absorption bands. The spectral channels of the O₂ A-band (759–770 nm) are strongly influenced by atmospheric absorption; therefore, this range of spectral channels was excluded from the SS calculation.

• Linear angular correction of reflectance. This correction aimed to normalize the solar illumination and viewing geometry of the observations across all PICS to a common reference geometry. The correction was applied independently to each spectral channel to account for the sensitivity of reflectance to variations in SZA and VZA. By doing so, the angular-induced variability in TOA reflectance was reduced according to Eq. (10).

  $$\begin{array}{}\text{(10)}& \begin{array}{rl}{\displaystyle }{R}_{\mathcal{S},\mathit{\lambda }}^{S_k}(t_n{)}_{\text{corr}}& {\displaystyle }=R_{\mathcal{S},\mathit{\lambda }}^{S_k}\left(t_n\right)-a_{\mathit{\lambda }}\left(\text{SZA}\right(t_n)-{\text{SZA}}_{\text{ref}})\\ {\displaystyle }& {\displaystyle }-b_{\mathit{\lambda }}\left(\text{VZA}\right(t_n)-{\text{VZA}}_{\text{ref}})\end{array}\end{array}$$

  where $R_{\mathcal{S},\mathit{\lambda }}^{S_k}(t_n{)}_{\text{corr}}$ denotes the corrected TOA reflectance, $R_{\mathcal{S},\mathit{\lambda }}^{S_k}\left(t_n\right)$ is the original TOA reflectance prior to correction, and a_λ and b_λ represent the sensitivities (slopes) of TOA reflectance with respect to SZA and VZA, respectively, and SZA_ref and VZA_ref are the reference angles, set to 45 and 0°, respectively.

4.2.2 Parameterization of the stability score for each site

After retrieving the reflectance time series, heat maps were generated for the PICS. For each spectral channel, the corresponding parameters were derived and are presented in Table 3. These parameters were selected to estimate the SS of each PICS (see Fig. 2). The reasons for choosing these parameters are given in Table 3. All selected parameters shared the same stability indicator, where lower values corresponded to higher temporal stability.

Table 3 Summary of statistical parameters used for temporal stability assessment of the PICS and their physical interpretation in evaluating the stability of reflectance time series.

Figure 2 Methodological framework outlining the key steps: retrieval of reflectance time series, derivation of temporal metric parameters for each spectral channel, computation of SS per channel, band_wise averaging of SS of spectral channels, and derivation of the final SS for each PICS site. The lower left block of the diagram shows the workflow of the SS computation for one spectral channel (spectral channel index 0). The process was iterative for all spectral channels. Later, the SS for bands and PICS were computed.

The selected parameters for each time series of a wavelength were normalized based on Min-Max scaling (Han et al., 2022), shown in Eq. (11). This normalization aims to prevent large numerical ranges from dominating the magnitude of SS. The scaled parameters were in the range of 0 to 1.

$$\begin{array}{}\text{(11)}& \stackrel{\mathrm{̃}}{F}={\displaystyle \frac{F-F_{\text{min}}}{F_{\text{max}}-F_{\text{min}}}}\end{array}$$

where F refers to the statistical parameters in Table 3. For each spectral channel, $\stackrel{\mathrm{̃}}{F}$ is the normalized parameter.

SS was computed based on Eq. (12), in which the weight of parameters was equal for all.

$$\begin{array}{}\text{(12)}& {\text{SS}}_{\mathcal{S},\mathit{\lambda }}^{S_k}=\sum _{k=\mathrm{1}}^K{w}_k{\stackrel{\mathrm{̃}}{F}}_k,\text{with}\sum _{k=\mathrm{1}}^K{w}_k=\mathrm{1},\text{and}{w}_k={\displaystyle \frac{\mathrm{1}}{k}}\end{array}$$

where K is the number of features.

The overall SS of a PICS was represented by the average of ${\text{SS}}_{\mathcal{S},\mathit{\lambda }}^{S_k}$ of the whole spectrum (see the illustrative diagram in Fig. 2).

$$\begin{array}{}\text{(13)}& \stackrel{\mathrm{‾}}{{\text{SS}}_{\mathcal{S}}^{S_k}}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{j=\mathrm{1}}^N{\text{SS}}_{\mathcal{S},\mathit{\lambda }}^{S_k}\end{array}$$

where j is the number of spectral channels, ${\text{SS}}_{\mathcal{S},\mathit{\lambda }}^{S_k}$ is the stability score of a spectral channel of a wavelength λ.

5 Results

5.1 GOME and SCIAMACHY reflectance time series of PICS

The reflectance time series of the PICS were retrieved under conditions of CF less than 0.25. Libya 4, one of the PICS endorsed by CEOS, WGCV, and IVOS, has been extensively used in previous studies for long-term stability assessments and sensor calibration (Lacherade et al., 2013; Barsi et al., 2018; Khadka et al., 2021; Helder et al., 2013). As a representative example, Libya 4 was selected for detailed visualization of the reflectance time series obtained from the GOME and SCIAMACHY sensors.

Figures 3 and 4 present the corresponding heatmaps of reflectance measurements from GOME and SCIAMACHY, respectively. At certain points in time for GOME (Fig. 3), higher values of the reflectance were observed across the whole spectrum (horizontal bright lines). In the UV, the value of the reflectance exceeded the average value due to instrument degradation, particularly towards the end of the GOME mission in 2011. The VIS and NIR spectral ranges were more uniform. In the NIR spectral range, O₂ A-band was visible with low reflectances around the center of the band. Similar behavior was observed for SCIAMACHY as shown in Fig. 4. The main difference was the absence of larger variations over long periods in the UV, suggesting better degradation correction compared with GOME.

Figure 3 Heatmap of the GOME reflectance observations (N=716) over Libya4. Three representative wavelengths from the UV, VIS, and NIR regions are highlighted with white dashed lines, and their time series are shown in Fig. 5.

Figure 4 Same as in Fig. 3 but for SCIAMACHY reflectance observations (N=568) and the time series of the selected wavelengths are shown in Fig. 6.

To gain further insight into the variability of the reflectance time series observed in the heatmaps, three wavelengths were selected. These wavelengths, located approximately at 330 (UV), 450 (VIS), and 770 nm (NIR), were indicated by white dashed lines in the heatmaps. Figures 5 and 6 show the corresponding reflectance time series at these wavelengths, derived from GOME and SCIAMACHY data, respectively.

Figure 5 The reflectance time series of the wavelengths: 329.96 nm in UV (blue), 449.98 nm in VIS (green), and 770.10 nm in NIR (orange) from GOME with (N=716).

Figure 6 The reflectance time series of the wavelengths 329.98 nm in UV (blue), 449.94 nm in VIS (green), and 770.04 nm in NIR (orange) from SCIAMACHY with (N=568).

The reflectance time series at the UV wavelength of GOME at 329.5 nm exhibited significant fluctuations, with approximately 47 % of the observations within ±10 % of the average value. A noticeable increase in reflectance was observed in 2001 and toward the end of the GOME mission due to degradation effects (see Fig. 5). Additionally, the reflectance time series exhibited three clearly separated clusters, likely resulting from the dependence of reflectance on VZA. In contrast, the selected wavelengths in the VIS (436.47 nm) and NIR (736.19 nm) regions demonstrated greater consistency, with 87 % and 96 % of the observations, respectively, remaining within ±10 % of their average values. Among them, the VIS wavelength showed a particularly more consistent pattern.

For SCIAMACHY, the average reflectance values at the corresponding wavelengths were comparable to those from GOME (see Fig. 6). A similar fluctuation pattern was observed in the UV band (329.98 nm), about 58.3 % of observations were within ±10 %. Additionally, sinusoidal-like patterns were more apparent in the SCIAMACHY UV time series than in GOME. Unlike the GOME UV reflectance time series, the scatter distribution did not show distinct clustering. The reflectance time series for the selected VIS (449.94 nm) and NIR (770.04 nm) wavelengths appeared even more consistent than those of GOME, with no evident degradation effects over the mission period.

5.2 Temporal metrics from GOME and SCIAMACHY reflectance time series

The derived temporal metrics of the GOME and SCIAMACHY reflectance time series of each spectral channel from the investigated PICS are presented in Figs. 7 and 8. These results offered a comprehensive comparison of statistical features across the PICS sites and spectral channels.

Figure 7 Temporal metrics of GOME reflectance time series of the investigated PICS include mean, standard deviation, coefficient of variation, interquartile range, slope, skewness, and kurtosis. Each subplot consists of three boxes, which refer to the wavelength ranges in UV/VIS/NIR.

Figure 8 Same as in Fig. 7 but for SCIAMACHY reflectance time series.

The mean of the GOME reflectance time series in Fig. 7a demonstrated relatively low reflectance values in UV and VIS and higher reflectance values near the O₂ absorption window (O₂ A-band). The standard deviation in Fig. 7b and IQR in Fig. 7d showed higher UV variability for most sites. However, tighter spreads of the IQR were observed in VIS (see Fig. 7d). The CV in Fig. 7c showed notably higher values and a wide range across spectral channels and sites. The CV values were suddenly increased in the O₂ A-band. This was likely due to the low mean reflectance value in the absorption window. The slope values in Fig. 7e showed a notable pattern of degradation across most of the sites, particularly for the short wavelengths in UV. Smaller degradation patterns were observed in the VIS and NIR ranges. GOME reflectance time series demonstrated large magnitudes of skewness in Fig. 7f and kurtosis in Fig. 7g, suggesting higher variability and asymmetry across the spectral channels, especially in the O₂ A-band in NIR, which could be due to the influence of the clouds and aerosols (Min et al., 2014).

The mean reflectance values of SCIAMACHY were similar to GOME and showed similar patterns across all sites in Fig. 8a. There were relatively low reflectance values in UV and VIS and higher reflectance values in NIR near the O₂ A-band. Low reflectance values were observed within the O₂ A-band (759–770 nm), indicating strong atmospheric absorption. The standard deviation in Fig. 8b and the IQR in Fig. 8d provided insights into the temporal variability of reflectance. Higher values in the UV range suggested greater fluctuations and noise, while lower values were in the VIS and on both sides of the O₂ A-band in NIR. The CV in Fig. 8c revealed differing levels of stability across spectral bands and among the PICS sites. In the UV, CV tends to increase with wavelength, while in the VIS, it remains relatively consistent. High CV values in the O₂ A-band in the NIR. Overall, the results suggested greater temporal stability in the VIS and spectral channels adjacent to the O₂ A-band compared to the UV. CV values among the sites range between 0.05 (5 %) and 0.20 (20 %). The slope values in Fig. 8e, especially in the UV, showed slight changes in the spectral signature of the sites. The oscillations in the slope of the UV channel were possibly caused by residual etalon variations. The dip around 350 nm was likely due to stronger degradation of the Band 2 detector, which was a known characteristic of the instrument. In contrast, many PICS sites exhibited slope values close to zero in the VIS and NIR, suggesting negligible long-term trends and stable reflectance behavior in these regions. The skewness values in Fig. 8f varied among sites and bands but remained consistent, reflecting differing degrees of asymmetry in the reflectance distributions. The VIS region exhibited more symmetric reflectance distributions compared to the UV. In the NIR, the absorption feature caused noticeable disruptions in symmetry. Similarly, kurtosis values in Fig. 8g, especially in the UV and NIR, suggested the presence of occasional outliers or extreme values. The VIS region, in contrast, showed fewer such anomalies.

The comparative analysis of temporal metrics for GOME and SCIAMACHY reflectance revealed different sensor performances and varying patterns across both PICS sites and spectral bands. Before assessing the stability of PICS, it is essential to first identify the most stable sensor data to be used in the analysis.

5.3 GOME versus SCIAMACHY: a comparative analysis

The averages of the metric temporal parameters of each spectral band and site are shown in Fig. 9. The average values were computed by taking the mean of the parameters across the wavelength ranges defined for each band in Table 2. Key findings from these temporal metrics for both sensors are presented in Table 4. It presents the comparison of GOME and SCIAMACHY for each extracted parameter of each band.

Figure 9 Average of metric temporal parameters for the investigated PICS. For each PICS site and sensor, three bars represent the metric values in the UV, VIS, and NIR bands (blue, orange, and green, respectively). Solid bars denote SCIAMACHY results, while hatched bars denote GOME.

Table 4 Comparison of the average of metric temporal parameters for the investigated PICS.

GOME and SCIAMACHY shared several characteristics, particularly in their mean reflectance profiles. Both sensors exhibited similar reflectance patterns across the spectral channels (Fig. 9a). While both sensors displayed high temporal variability in the UV band, they had less variability in VIS and the lowest in NIR (see Fig. 9b and c).

However, GOME consistently displayed higher temporal variability and drift indicators than SCIAMACHY. Specifically, GOME's standard deviation, IQR, and CV were relatively higher (Fig. 9b–d); its slope values suggested greater sensor drift (Fig. 7e); skewness was more strongly positive across most sites (Fig. 9f); kurtosis revealed more frequent extreme reflectance events (Fig. 9g). Collectively, these findings recommend adopting SCIAMACHY as the primary reference sensor for PICS stability assessments.

5.4 Reflectance correction for sun illumination and satellite viewing angles

To illustrate the effects of SZA and VZA on reflectance, the mean reflectance across all spectral channels was analyzed as an example. Figure 10a and b present scatter plots of the mean reflectance of the whole spectrum as a function of SZA and VZA, respectively. The black dashed lines represent linear regression fits and illustrate the overall trends in the observations. Consistent positive slopes were observed for all sites with increasing SZA, indicating an increase in reflectance with SZA. In contrast, reflectance exhibited consistent negative slopes with increasing VZA.

Figure 10 Mean reflectance of the whole spectrum of the PICS observations with respect to SZA and VZA. The scatterplot of the mean reflectance versus (a) SZA, (b) VZA. The fitted linear regression equations for the samples are shown in the black dashed lines (c) the percentage of the changes in mean reflectance per 10° changes in SZA (blue bars) and VZA (orange bars).

Figure 10c shows the percentage change in reflectance corresponding to a 10° change in the sun-illumination and viewing geometry. The mean reflectance typically varies by 0.5 % per 10° to 2.57 % per 10° increase in SZA, with some sites, such as Arabia1 and Niger1, showing lower sensitivity of approximately 0.5 %. The effect of VZA was more pronounced, exhibiting an inverse relationship with reflectance. The corresponding absolute variations in reflectance ranged from 2 % per 10° to 6 % per 10° change in viewing angle.

The regression analysis revealed a systematic dependence of mean reflectance on both solar illumination and satellite viewing geometry. The magnitude of changes in reflectance was not the same across the PICS. Consequently, the reflectance observations over the PICS were corrected for the effects of both SZA and VZA. To account for spectral variability, the regression was applied separately to each spectral channel, ensuring that angular dependencies were corrected on a channel-by-channel basis.

As an illustrative example, the reflectance differences for the Libya4 site (corrected – original) are provided in the Supplement (see Fig. S1). The effect of the correction was strongest in the UV spectral channels, weaker in the VIS channels, and smallest in the NIR region. The correction applied to the reflectance time series exhibited a clear seasonal dependence, with the largest correction values for the reflectance time series of the summer months (May–September).

5.5 Stability score for PICS

The SS for the PICS sites were computed based on corrected SCIAMACHY reflectance observations. The SS was calculated for each spectral channel, excluding the O₂ A-band. The scores were normalized and visualized in Fig. 11 with lower values indicating greater stability.

Figure 11 The heatmap of the SS of each spectral channel of UV, VIS, and NIR spectral channels. The O₂ A-band is excluded from SS computation and appears as a white rectangle in the NIR panel.

Some sites demonstrated stable performance across nearly all spectral channels. With the most stable sites being Sudan1, Arabia2, and Libya4. In contrast, other sites showed differences when comparing the stability of the different spectral bands. For instance, Arabia3, Arabia2, Libya3, and Niger2 showed higher stability in the UV band compared to the VIS band. Conversely, sites such as Algeria1, Algeria2, Algeria3, Libya4, and Sudan1 demonstrated greater stability in the VIS band than in the UV. A few sites, including Mali1, Mauritania1, Mauritania2, and Niger1 displayed generally low stability across all spectral regions. Additionally, high SS values were observed across multiple sites, appearing as vertically homogeneous patterns. This effect was particularly pronounced in the UV spectral region, likely due to the strong sensitivity of UV bands to atmospheric influences, such as cloud contamination and aerosols.

Based on the band-wise ranking shown in Fig. S2, the SS values of the PICS span a wide range across the spectral domains. In the UV range, SS varies from 0.11 to 0.75, while in the VIS range, it lies between 0.25 and 0.70. The NIR range exhibits SS values between 0.16 and 0.78. On average, the SS values for the UV and VIS bands are 0.45 and 0.50, respectively, which are higher than the average value for the NIR band (0.40). This indicates that PICS stability exhibits a clear dependence on spectral band. As an illustrative example, the Libyan sites (Libya1, Libya2, Libya3, and Libya4) exhibit SS values of 0.50, 0.38, 0.36, and 0.33, respectively (see Fig. S3). The high SS of Libya1 can be attributed to its being relatively higher SS in the UV and VIS bands compared to the other Libyan sites, as the SS represents a combined measure across all spectral bands.

The influence of geographical proximity is linked to the high potential of similar surface characteristics, which is evident in the band-wise ranking of the PICS, particularly in the NIR band, excluding the O₂ A-band. Several geographically close sites showed similar SS values, including Algeria3 and Algeria4, Arabia1 and Arabia2, Mauretania1 and Mauretania2, Libya4 and Egypt1, and Libya3 and Libya2 (see Fig. S2). These geographically driven groupings – most clearly observed in the NIR band – suggest that sites within similar regional and geomorphological settings tend to exhibit comparable temporal characteristics. In contrast, this clustering effect was less pronounced in the VIS and UV bands.

6 Discussion

6.1 Framework

To ensure the robustness and comprehensive assessment of the reflectance time series and the comparison between sensors, additional statistical descriptors were included. Besides the basic statistical parameters and CV, other parameters like IQR, slope, skewness, and kurtosis were included in the assessment. Skewness, kurtosis, and IQR parameters are used to identify the level of symmetry in the observations, reflect the presence of outliers, and quantify the spread of the central portion of the data. The slope serves as an indicator of the degradation in the sensor performance.

With only basic statistical parameters (mean, standard deviation, and CV), there were no significant differences observed for the spectrometer samples. When additional parameters such as slope, skewness, and kurtosis were considered, meaningful differences became apparent (see Fig. 9). These high-order statistical features provided greater sensitivity to distributional characteristics, subtle anomalies, and temporal trends, thereby offering a more robust basis for distinguishing the performance and stability of the spectrometers. This multi-parameter approach contributed to a more robust and resilient framework for sensor performance evaluation, radiometer drift monitoring, and identifying the constrained and reference sensors for cross-calibration.

6.2 Reflectance time series analysis

The derived indicators from the reflectance time series of GOME and SCIAMACHY revealed different statistical patterns among the PICS and across spectral channels. These differing patterns indicated that the initial hypothesis – that atmospheric conditions across the sites are relatively similar and that all PICS should exhibit similar behavior – may not hold in all cases. Therefore, PICS should be regularly investigated (Tuli et al., 2019). This finding underscored the need to quantify and compare the stability of each PICS site individually and across different spectral bands. The metric parameters derived from reflectance observations from both sensors consistently showed that deviation, variation, asymmetry, and degradation were highest in the UV bands, lower in the VIS bands, and lowest in the NIR bands. These results were expected for the following reasons: (i) TOA reflectance was used, which includes contributions from both surface and atmospheric reflectance; (ii) the UV bands are highly sensitive to atmospheric processes and aerosols (Chatzopoulou et al., 2025; Li et al., 2012); (iii) the VIS bands are influenced by Rayleigh scattering and aerosol variability (He et al., 2021); and (iv) the NIR bands (excluding the O₂ A-band) are less affected by atmospheric processes and are more sensitive to surface reflectance; therefore, it serves as a baseline for atmospheric correction algorithms (Mishra et al., 2020).

6.3 Comparative analysis of GOME versus SCIAMACHY

A comparative analysis of key statistical parameters revealed that SCIAMACHY exhibited more statistically uniform and stable reflectance measurements than GOME. These indicators suggested lower temporal variability and fewer extreme values in SCIAMACHY data, making it a more robust choice. Noël et al. (2007) investigated the degradation across the spectral channels of SCIAMACHY and found that even after reaching the nominal lifetime of the mission in 2007, the SCIAMACHY instrument was still in good condition. However, the same study demonstrated an increased degradation in UV channels.

The discrepancy in the derived statistical parameters between GOME and SCIAMACHY might be partly attributed to sensor degradation effects. Tilstra and Stammes (2006); Coldewey-Egbers et al. (2008) referred to the severe degradation of the GOME instrument in the UV due to scan mirror degradation. Overall, SCIAMACHY's enhanced design and calibration capabilities resulted in better long-term radiometric stability compared to GOME. Hilbig et al. (2019) and Krijger et al. (2014) presented a degradation correction approach based on SCIAMACHY's internal white light measurements in combination with direct solar measurements. Consequently, SCIAMACHY was adopted as the reference sensor for evaluating the temporal stability of the PICS. Using SCIAMACHY observations in this paper ensured that the derived SS were based on a more stable and consistent dataset, enhancing the reliability of the site rankings and subsequent cross-calibration of both sensors by Owda et al. (2026).

6.4 Refinement reflectance observations for reliable estimation of the stability scores

The reflectance observations collected over PICS exhibited varying degrees of dependence on the SZA and VZA across the sites. An increase in SZA results in a longer atmospheric path length, leading to enhanced scattering and, consequently, TOA reflectance. This dependence on SZA is consistent with previous studies; for example, Wang et al. (2022) reported that reflectance tends to increase with both atmospheric aerosol optical depth and SZA. These findings indicated that a significant portion of the observed reflectance variability could be attributed to the sensitivity of TOA reflectance to changes in SZA. Differently, the reflectance time series as a function of VZA exhibited an inverse relationship, which may be attributed to degradation of the scan mirror angle over time, particularly in the UV band.

The correction applied to account for reflectance dependence on variations in SZA and VZA exhibited a clear seasonal pattern, with maximum corrections occurring during summer. Using a reference geometry with SZA=45°, the typically lower SZA values during summer led to increased reflectance after correction to the reference geometry compared with winter. The magnitude of the correction was wavelength dependent, with the largest corrections observed in the UV bands. This behavior was likely associated with enhanced atmospheric transport of aerosols and dust during summer. Previous studies reported a pronounced seasonal dust cycle, characterized by higher dust concentrations in the summer months (Rezaei-Nokandeh et al., 2025). Both the central Sahara Desert and the Arabian Peninsula experienced their highest aerosol optical depth during summer (Logothetis et al., 2021). Accounting for these effects was therefore essential for modeling SS and for ensuring comparable values across sites, without the influence of differing illumination and viewing geometries.

6.5 Stability of PICS

SS served as an indicator of temporal variability across the spectral channels. The SS for PICS were estimated based on groups of parameters that directly quantified variability, distribution, and temporal trends. The features describing temporal variability in reflectance were normalized to ensure comparability across sites. The SS was primarily determined by band characteristics. Overall, the values of SS were lower (more stable) in the NIR than in the VIS and UV bands. At longer wavelengths, contributions from surface properties – such as mineral composition, soil texture, surface roughness, and moisture – become increasingly significant (Fang et al., 2018), increasing the sensitivity of the metric to surface homogeneity and underlying geophysical conditions.

The band-wise ranking of the PICS showed strong clustering based on geographical proximity, which can be attributed to similar soil mineralogy and morphological conditions. In contrast, less clustering by geographic proximity was observed in the UV and VIS bands. This likely reflects regional differences in meteorological and atmospheric conditions, as these bands are more strongly influenced by atmospheric effects.

Overall, the results demonstrated a wide range of SS across spectral bands and sites. Our findings confirmed that Libya4 – one of the recommended sites for sensor cross-calibration – was among the most stable PICS. Some sites, such as Sudan1 and Arabia2, exhibited lower scores (more stable) than Libya4. These three sites are, therefore, recommended for further cross-calibration activities.

Differing SS across the PICS could refer to the following reasons:

• Surface characteristics. Site properties affect radiometric stability; PICS are not fully invariant. Relevant factors include sand topography, grain size (Niro et al., 2021), composition, color, material type, and dune displacement (Govaerts, 2015). The sun's azimuth relative to dune alignment influences reflectance, with mean values depending on the area size (Govaerts, 2015).

• Meteorological conditions. Strong winds can transport dust over large areas. Saharan dust storms in June 2020 caused high aerosol optical depths to extend to the USA within 4 d (Francis et al., 2020). Such events particularly affect short wavelengths (blue/visible) due to aerosol scattering.

• Atmospheric conditions. The derived metrics used to compute the SS were based on TOA reflectance. Consequently, the SS remain influenced by atmospheric conditions, particularly in the UV and VIS spectral bands, where atmospheric scattering and absorption effects are more pronounced. Contributions from aerosols, molecular (Rayleigh) scattering, and trace gases can introduce additional variability in the TOA signal, which may not be fully separable from surface-related changes. As a result, the computed SS reflect the combined effects of surface properties and atmospheric variability rather than purely surface stability. Nevertheless, the analysis focuses on relative temporal behavior and inter-site comparisons, the impact of atmospheric effects is assumed to be comparable across the selected PICS, allowing meaningful assessment of their relative stability.

7 Conclusion

This paper presented a statistical framework based on multiple indicators to compare different sensors, monitor sensor performance, identify sensor-related issues, determine reference and constrained sensors for cross-calibration, and investigate the temporal variability of spectral bands and PICS based on TOA reflectance. The framework incorporated a wide range of parameters, including traditional time series metrics such as the mean, standard deviation, and CV. In addition, it included additional indicators that provide insights into the distribution shape of observations, the presence of anomalies, and temporal trends.

The GOME and SCIAMACHY spectrometers were used in this study due to their similar sensor characteristics and mission objectives, as well as overlapping spectral coverage. Additionally, their missions provided approximately 10 years of overlapping observation data, with a temporal offset of about 30 min in their overpass times. A decade of reflectance measurements from both sensors was retrieved under nearly clear-sky conditions over PICS.

The results highlighted instrumental issues in GOME, including significant degradation observed in 2001 and irregular observation patterns toward the end of its operational life. These findings agreed with previous studies that reported pronounced degradation of GOME, particularly in the UV spectral band. Moreover, statistical indicators revealed high positive values of skewness and kurtosis, indicating an asymmetrical distribution and the presence of numerous outliers. In contrast, SCIAMACHY observations were found to be more uniform and stationary, with substantially lower trends compared to GOME. As a result, SCIAMACHY was considered the reference sensor in this study for assessing the stability of PICS.

In terms of site stability, areas larger than the spatial extension of PICS were considered due to the large footprint of the spectrometers. Furthermore, this study aimed to establish a foundation for spectrometer cross-calibration, for which the temporal stability of areas larger than the PICS is a critical factor. The analysis revealed considerable variation among the PICS. While some sites exhibited high temporal stability across all spectral channels, others showed inconsistent behavior across bands, and a subset of sites displayed significantly lower stability, up to two to three times less stable than the most consistent sites. Additionally, the SS exhibited a clear band-wise dependency, with higher stability observed in the NIR band (excluding the O₂ A-band) and lower stability in the VIS and UV bands. These findings emphasize the importance of selecting suitable PICS for radiometric calibration, taking into account the spectral band of interest.

The identification of the most stable PICS has broader implications for the satellite remote sensing community. Highly stable PICS can serve as reliable reference targets for monitoring long-term sensor performance, assessing radiometric consistency, and supporting cross-calibration activities. Improved radiometric stability at the calibration level directly propagates to higher-level satellite products. In particular, the retrieval of atmospheric trace gases depends critically on the accurate spectral structure of TOA reflectance. For example, Differential Optical Absorption Spectroscopy (DOAS) retrievals typically require a relative point-to-point accuracy on the order of 10⁻⁴. This stringent requirement highlights the importance of carefully selecting stable reference sites for tracking sensor performance, radiometric calibration, and cross-calibration. Furthermore, periodic evaluation of these sites is crucial to ensure their long-term stability. The follow-up paper performs the cross-calibration of GOME and SCIAMACHY based on the most stable sites of this research paper (Owda et al., 2026).

The findings of this study indicate that the initial hypothesis does not hold uniformly across all sites, as atmospheric conditions can vary significantly between desert regions. We recommend that future research focus on separating the atmospheric reflectance component from the TOA reflectance. Comparative statistical analyses of the NIR relative to the UV and VIS spectral ranges can provide insight into the stronger atmospheric influence in the UV and VIS. Under this framework, NIR observations may serve as a reference to correct UV and VIS reflectance. In addition, incorporating auxiliary parameters – such as aerosol optical depth and trace gas column densities derived from the corresponding Level-2 products – into the statistical analysis could further improve the robustness of the methodology.