Spectral decomposition of remote sensing reflectance variance due to the spatial variability from ocean color and high-resolution satellite sensors

Abstract. The variability of the remote sensing reflectance, Rrs, now routinely retrieved from ocean color (OC) and high spatial resolution sensors, is often used to characterize water variability due to changes in inherent optical properties of the water body. At the same time, Rrs is partially variable because of uncertainties in its retrieval in the process of atmospheric correction. Using data from SNPP-VIIRS and Landsat-8 OLI sensors, the contribution of the main components to the variance of Rrs due to its spatial variability is determined based on a model in which variances were considered proportional to the mean values of the corresponding components. It is shown that there is practically no spatial variability in the open ocean waters and water variability is proportional to the spatial resolution of the sensor in coastal waters. Variances due to surface effects, inaccuracies of aerosol models, and sunglint can contribute significantly to Rrs variance, which characterizes Rrs spatial variability, with variances due to the water variability itself often being significantly smaller.


Introduction
Ocean color (OC) is indicative of ocean health and biochemistry and for that reason is listed as an essential climate variable (ECV). 1 The color of a water body is determined by scattering and absorption by pure water and its natural constituents, such as phytoplankton, other particles, and colored dissolved organic matter (CDOM). 2 Some phytoplankton species form harmful algal blooms, which can negatively affect human and marine life and often have severe repercussions on a range of industries. 3 Typically, <10% of the top of the atmosphere (TOA) radiance is due to the water signal at sea level, 4 with the remainder originating from scattering processes in the atmosphere and reflections of Sun and sky on the wave-roughened water surface. It is paramount to accurately estimate radiances at the surface level from the ones at the TOA as uncertainties propagate into the retrieval of water parameters, characteristics of in-water particulates, concentrations of chlorophyll-a, and detection of algal blooms. 5 Atmospheric correction uncertainties stem at least partially from the estimation of aerosol models 6,7 and air-water interface corrections, [8][9][10] with the latter due to sky and Sun light reflections at the wind-roughened air-water interface. Glint correction uncertainties often have a stronger impact on retrievals in coastal waters with low water reflectance values in the blue bands. [11][12][13][14] Currently, high uncertainties in blue reflectance observations are widely acknowledged. 5,14,15 However, little is known about the specific dependencies of uncertainties *Address all correspondence to Alexander Gilerson, gilerson@ccny.cuny.edu concerning their spectral and scaling behavior in various water areas and their dependence on meteorological conditions. Uncertainties of the remote sensing reflectance, R rs , are estimated by comparison of the water parameters determined from the satellite imagery with the "true" values, which can be determined, for an example, in very uniform clear waters in which all water parameters can be connected to the concentration of chlorophyll-a, [Chl]. 16 A second approach is to compare data from satellite with field measurements from the towers in the ocean, the Aerosol Robotic Network for Ocean Color (AERONET-OC 17,18 ) sites or from ships. 19 Such comparisons can be carried out in a wide range of waters; however, they are associated with multiple additional uncertainties, which are related to the quality of field measurements themselves, 20 water variability inside pixels, and time difference between in situ and satellite data. In another approach that uses Monte Carlo (MC) simulations for sea-viewing wide field-of-view sensor (SeaWiFS) observations, 21 the retrieval process for the remote sensing reflectance, R rs , was repeated 1000 times, and uncertainties in R rs were then estimated as the "standard deviation of the 1000 perturbed R rs retrievals in each band." This derived uncertainty was interpreted "as the precision of the R rs retrieval due to instrument noise." The standard deviation of R rs due to its spatial variability in the specific area, σ spat , represents the changes related to the variability of water inherent optical properties and together with the coefficient of variation (CV) is of interest to the OC community, especially as a function of the satellite spatial resolution or ground sampling distance (GSD). 22 However, σ spat is contaminated by other components of R rs uncertainites and should be rectified before being considered a proxy for water variability. σ spat is also a part of the full R rs uncertainty and can be considered the precision of the R rs retrieval due to spatial variability for the specific type of the atmospheric correction process.
In this work, we focus on the spectral behavior of σ spat and its components and their dependence on GSD through the analysis of satellite data from the Visible Infrared Imaging Radiometer Suite (VIIRS) on the Suomi National Polar-orbiting Partnership (SNPP) platform and Landsat-8 Operational Land Imager. First attempts in this direction were made through the analysis of VIIRS satellite data, 23 demonstrating dependence of spectra and magnitude of σ spat on GSD, which were attributed to the surface and water conditions. This work will expand on that study by considering a broad set of parameters involved in the atmospheric correction process.

Theoretical Considerations
The main radiometric quantity in the processing of satellite data is remote sensing reflectance, R rs , which is defined as the ratio of the water-leaving radiance to the downwelling irradiance at the sea surface, R rs ðλÞ ¼ L w ðλÞ∕E d ðλÞ, where L w ðλÞ is the water leaving radiance, E d ðλÞ is the downwelling irradiance, and λ is the wavelength. At the TOA, the total radiance, L t ðλÞ, is represented as 6 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 2 8 3 L t ðλÞ ¼ L r ðλÞ þ L a ðλÞ þ L g ðλÞ þ tðλÞL wo ðλÞ; (1) where L r ðλÞ is the radiance due the Rayleigh scattering including surface effects, L a ðλÞ is the radiance due to the aerosol scattering, L g ðλÞ is the radiance due to the sun glint from the water surface, L wo ðλÞ is the water leaving radiance in the process of detection, and tðλÞ is the diffuse transmittance of light from the water surface to the TOA. L t ðλÞ has uncertainties due to all of these components and sensor noise. In the process of the retrieval of the water leaving radiance E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 1 9 2 L w ðλÞ ¼ ðL t ðλÞ − L r ðλÞ − L a ðλÞ − L g ðλÞÞ∕tðλÞ: In addition, L r ðλÞ can be divided into the radiance from the Rayleigh scattering in the atmosphere, L R ðλÞ, and the radiance due to the reflectance from the ocean surface, L surf ðλÞ E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 1 3 7 L r ðλÞ ¼ L R ðλÞ þ tðλÞL surf ðλÞ; (3) where L surf ðλÞ ¼ L sky ðλÞ Ã ρ, L sky ðλÞ is the sky radiance and ρ is the reflectance coefficient from the water surface. In the satellite atmospheric correction procedure, averaged surface effects are considered in the vector radiative transfer (VRT) equations, which are based on Cox-Munk distributions. 8,24 However, each satellite image captures a specific snapshot of the ocean, in which the actual spatial average of the light field reflected from the wave facets may not exactly match the average predicted by the VRT model. The actual signal may have its own features due to instantaneous water and atmospheric conditions, spatial scales in the area, or simplifying assumptions made within the VRT model, as Cox-Munk model is not necessarily valid for waters in coastal areas. Uncertainties from the abovementioned components need to be taken into account. So, normalizing by the downwelling irradiance, E d ðλÞ, uncertainty in remote sensing reflectance σ in sr −1 is determined from E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 6 3 9 where σ 2 t , σ 2 R , σ 2 a , σ 2 g , σ 2 surf are variances of normalized L t ðλÞ, L R ðλÞ, L a ðλÞ, L g ðλÞ, L surf ðλÞ and σ 2 water ðλÞ and σ 2 noise ðλÞ are variances due to water variability and sensor noise respectively. Variances for the quantities at TOA σ 2 t , σ 2 R , σ 2 a , σ 2 g are divided by t 2 in accordance with Eqs. (1) and (2); σ 2 noise has 1∕t 2 in its definition. 25 We assumed that all standard deviations (except σ noise ) in Eq. (4) as a first approximation, are spectrally proportional to the corresponding mean values of normalized radiances with k as proportionality coefficients (k t , k R , k a , k g ,k s , k R rs , and k n are the fitting coefficients for L t ðλÞ, L R ðλÞ, L a ðλÞ, L g ðλÞ, L surf ðλÞ normalized by E d ðλÞ, R rs ðλÞ and σ noise respectively): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 5 0 5 and in Eq. (6a), F 0 ðλÞ is the extraterrestrial irradiance, τ R ðλÞ is the Rayleigh optical thickness, θ is the sensor zenith angle, and Θ is the scattering angle, which is the angle between solar and viewing directions; in Eq. (6b, τ a ðλÞ is the aerosol optical thickness, P a is the scattering function for aerosols assumed as P a ¼ 0.3 (without considering its possible spectral dependence) based on the analysis for two AERONET-OC stations; 26 in Eq. (6c), T 0 ðλÞ and TðλÞ are the direct transmittance coefficients for TOA to surface and surface to TOA, respectively, and 0.005 is the threshold for glint detection 9 , L GN , in sr −1 ; in Eq. 6(d), θ 0 and t 0 ðλÞ are the Sun zenith angle and the corresponding diffuse transmittance from the TOA to the water surface, respectively; in Eq. (6e), a representative normalized sky reflectance, S ¼ L sky ∕E d , was simulated by the VRT RayXP code 27 for the Sun zenith angle θ 0 ¼ 42 deg at a viewing zenith angle of 40 deg, aerosol optical thickness at 443 nm AOT 443 nm ¼ 0.13 and Angstrom coefficient 1.0, parameters that are very close to the average ones for all observations and are used throughout this study. Sky radiance, L sky , differs from scene to scene, whereas the normalized radiance, L sky ∕E d , has much less variability; the reflectance coefficient was considered to be ρ ¼ 0.025, which is the typical reflectance coefficient at a 40 deg viewing angle. It should be emphasized that the spectral shapes of the main components in Eqs. (6a)-(6e) are of primary interest; changes of values that were assumed constant in the model do not affect these shapes or the contribution of the corresponding uncertainties from these components to the total σ 2 ðλÞ. The glint and noise components are included in the expression for the total variance [Eq. (4)] but were not included in the total radiance L t ðλÞ in Eq. (6f) since ideally their values should be zero.
The model of Eqs. (5) and (6) can be used to analyze contributions of spectral components to the total σ 2 ðλÞ if the latter can be determined from the comparison of satellitein situ data. Our current work is primarily focused on the contributions of the same spectral components but only to the variance σ 2 spat referring to the R rs ðλÞ spatial variability in the area, which can characterize water variability, so σðλÞ was replaced by σ spat ðλÞ in Eq. (5) and in the further analysis.
σ noise was determined based on typical radiance at TOA L t ðλÞ, signal-to-noise ratio (SNR), and transmittance coefficients t 0 ðλÞ for the TOA to the surface and tðλÞ for the surface to TOA. 25 σ noise for VIIRS, SeaWiFS, Moderate Resolution Imaging Spectroradiometer (MODIS), and Landsat 8 Operational Land Imager (OLI) with data taken from Refs. 16, 28, and 29 are shown in Fig. 1. A preliminary study 23 showed that in moderately clear waters σ spat ðλÞ was close to the standard deviation spectrum, σ F , determined by Franz and Karaköylü 21 for SeaWiFS, and this spectrum was added for comparison.
Thus, in terms of the single pixel processing, noise can be a significant contributor to σ spat ðλÞ, but with the often-used representation of data averaged over multiple pixels, this contribution is reduced inversely proportional to the square root of the number of pixels.
Mean spectra of R rs ðλÞ and σ spat ðλÞ, quantities from Eq. (6) were determined by processing of VIIRS data in the specific area with the wind speed W in each of the ranges W < 3 m∕s, 3 < W < 5 m∕s, and W > 5 m∕s. With these parameters calculated for a scene, a nonlinear least-squares fit was carried out in MATLAB using the Levenberg-Marquardt technique 30,31 to determine the respective values of k coefficients for each scene or set of scenes, taking advantage of the differences in the spectra of σ spat ðλÞ components in Eq. (5). These coefficients are assumed to be representative of the respective contributions of the atmospheric, water, and sky reflection components to the observed spatial scale variance σ 2 spat ðλÞ of the R rs since R rs is the result of the processing through the atmospheric correction. By this operation, the actual contribution of σ water ðλÞ to σ spat ðλÞ is determined at different spatial resolutions since σ water ðλÞ better characterizes water optical variability than σ spat ðλÞ. Since, as was mentioned before, in the clear waters, the σ spat ðλÞ spectrum was very close to σ F ðλÞ, 21 the latter will be shown as a benchmark for the results given below that are related to σ spat ðλÞ. for three AERONET-OC site areas: the Long Island Sound Coastal Observatory (LISCO), WaveCIS site in the Gulf of Mexico, University of South California (USC) site in California waters, and the area of Marine Optical Buoy (MOBY) in Hawaii. NASA Level 2 data files for VIIRS include geophysical products of the atmosphere and ocean, such as aerosol optical thickness (AOT), remote sensing reflectance, R rs ðλÞ, and the level 2 quality flags. The Sun zenith angle, sensor viewing angle, and sensor azimuth angle were used to determine radiances and irradiance in Eq. (6). The VIIRS data have a nadir resolution of 750 m. Pixels used for matchup comparison were averaged over three spatial resolutions (3 × 3, 5 × 5, and 7 × 7 pixel boxes), centered at the in-situ observation location, 35 and the standard deviation between pixels was recorded. Pixels flagged by at least one of the following conditions were excluded: land, cloud, failure in atmospheric correction, stray light (except for LISCO), bad navigation quality, high or moderate glint, negative Rayleigh-corrected radiance, negative water-leaving radiance, viewing angle greater than 60 deg, and solar zenith angle greater than 70 deg.

Landsat-8 OLI Data
The Landsat-8 OLI has a resolution of 30 m and a field-of-view of 15 deg (AE7.5 deg from nadir). To analyze the dependence of uncertainties on various scales, the aquatic reflectance product was used. OLI bands are centered at 440, 480, 560, and 660 nm wavelengths. The L2 data were downloaded from the USGS Earth Resources Observation and Science Center Science Processing Architecture on Demand Interface website 36 for the period 2013 to 2019. R rs ðλÞ was calculated by dividing the dimensionless aquatic reflectance by π and was taken for those imageries that pass the selection by the level 2 quality flags. Pixels flagged by at least one of the following conditions were excluded: land, cloud, failure in atmospheric correction, bad navigation quality, high or moderate glint, sea ice, viewing angle greater than 60 deg, and solar zenith angle greater than 70 deg. For each spatial resolution, R rs ðλÞ was averaged, and the standard deviation was calculated. Nine spatial resolutions were implemented for each imagery with the following square pixel box scale: 3, 7, 19, 31, 35, 51, 75, 125, and 175 pixels (90, 210, 570,  930, 1050, 1530, 2250, 3750, and 5250 m, at nadir). Atmospheric correction for OLI 37,38 was implemented across 7 × 7 pixels to minimize noise effects, which need to be considered in the dependence of various parameters on the spatial resolution.
In both VIIRS and OLI processing, R rs ðλÞ spectra with R rs ð412Þ > 0.006 water were considered open ocean; otherwise, they were considered coastal water area cases. Spectra with R rs ð412Þ < 0, which are typically due to the very inaccurate aerosol model in the atmospheric correction, were not included in the processing. Such cases require advanced atmospheric correction procedures, 39 a discussion of which is outside the scope of this work. The areas of study and scenes analyzed from OLI imagery are shown in Fig. 2.

VIIRS Data Analysis
There are seven k coefficients in Eq. (5) that need to be determined from the fitting procedures, and VIIRS has only five bands in the visible part of the spectrum. Four terms were considered to be main potential contributors to σ spat ðλÞ: σ R ðλÞ, σ a ðλÞ, σ surf ðλÞ, and σ water ðλÞ; it was assumed that σ t ðλÞ was represented by its components and that the impact of two other terms, σ g ðλÞ and σ noise ðλÞ, was analyzed by adding one of them as a fifth term. It was found that the impact of noise σ noise ðλÞ was small and the contribution of the glint component σ g ðλÞ was much more pronounced. Further results are shown for the fitting with σ g ðλÞ being the fifth component. Typical spectra for all components involved normalized to the value at 412 nm are shown for the open ocean and coastal water station in Fig. 3. It should be emphasized that in Eq. (4) σ 2 R and σ 2 a are divided by the transmittance coefficient and σ 2 surf is not, so a combination of Rayleigh and aerosol scattering contributions has a normalized spectrum that is different from the spectrum of surface effects, which is proportional to the sky radiance. For the open ocean, the MOBY site surface effects spectrum in Fig. 3(a) is very close to the normalized spectrum of σ spat ðλÞ and, as shown in Fig. 4, is practically the only component represented in the fitting procedure, with the σ surf ðλÞ and σ fit ðλÞ being undistinguishable for all three wind speed ranges and σ fit ðλÞ being the best fit of spectral components into σ spat ðλÞ with k coefficients. These spectra are very close to the spectrum σ F ðλÞ. 21 The surface component is also dominant at the USC site open ocean waters with a small contribution of σ R ðλÞ and σ a ðλÞ components. The corresponding spectra of R rs ðλÞ in For the coastal waters in Fig. 5, there are some contributions from the Rayleigh and aerosol components (USC and WaveCIS sites), but the main effects are from the surface sky reflection, glint, and water variability. The changes of the k coefficients related to all of these effects with the spatial resolution will be shown below together with Landsat data.

Landsat Data Analysis and Comparison with VIIRS Data
NASA atmospheric correction for Landsat OLI is similar to the atmospheric correction for VIIRS, 37,38 and OLI fine spatial resolution permits an expanded analysis of σ spat ðλÞ dependence on GSD in a broad range. However, OLI has almost vertical viewing, whereas VIIRS viewing angles are in the range of 0 deg to 56 deg, which can impact Sun and sky glint contributions. In terms of the application of Eqs. (4)-(6) of the proposed model to σ spat ðλÞ from OLI, the main difficulty is that OLI has only four bands (whereas VIIRS has five bands), which limits simultaneous analysis of the critical parameters' contributions. After preliminary analysis, σ R ðλÞ was excluded from the consideration, and only contributions of aerosol, σ a ðλÞ, surface effects, σ surf ðλÞ, sun glint, σ g ðλÞ, and water variability, σ water ðλÞ, were analyzed. Examples of the fitting procedures for open ocean waters near the Hawaii site and WaveCIS coastal site for three spatial resolutions are shown in Figs. 6 and 7. σ spat ðλÞ is small at low resolutions and gradually increases for higher resolutions being very close to σ F ðλÞ, with main contributions being from surface effects and aerosol components in the open ocean waters.
At the WaveCIS site, there are clear contributions from all components, with increasing effects due to water variability proportional to the R rs ðλÞ spectrum for higher GSD. Mean values of normalized radiances for different components in Eq. (6) vary significantly, so k coefficients are given in Figs. 6 and 7 to demonstrate the contribution of components associated with these coefficients. Since there were no clear differences in the dependence of k coefficients on wind speed, changes of the fitting k coefficients as a function of GSD averaged over all wind speeds are shown in Fig. 8 for OLI.
The main outcomes from the application of Eqs. (5) and (6) of the model to OLI data at all considered sites are as seen in Figs. 8 and 9. In the open ocean waters (Hawaii area), the main contributions come from the surface effects gradually increasing from small values at small GSDs to k s ¼ 0.1 to 0.12 after about 1500 m. Recalling that the surface effect itself was considered the sky normalized radiance multiplied by the average surface reflectance coefficient  0.025, k s ¼ 0.1 corresponds to 10% changes in this value. There were also contributions from the aerosol component through k a ¼ 0.008 to 0.015, with its relative impact being visible in Fig. 6. There is not any spatial variability in Hawaii waters until GSD of about 1500 m, and it is very small at larger distances in terms of the k R rs coefficient. Glint effects were not visible in open ocean waters. In coastal waters, surface effects are smaller than in the open ocean and vary from site to site with average k s being about 0.06 and larger contributions typically being from the aerosol component k a ≈ 0.02, with some exceptions. As a reminder, in Eq. (6c), the glint mean value was determined based on L GN ¼ 0.005, so k g ¼ 0.02 to 0.05 represents 2% to 5% of this value. The most prominent effect is seen from the water variability with k R rs changing approximately proportionally to GSD and reaching about 7% to 14% at about 5000 m.
Due to different NIR wavelengths being used in the atmospheric correction models for OLI and VIIRS and different viewing angles, a good match among the k s , k a , k g coefficients for these two sensors was not expected. To analyze the consistency of σ spat components for two sensors, a combined variance σ 2 sag ("s" for the surface, "a" for aerosols, and "g" for glint in the subscript) was calculated as σ 2 sag ¼ σ 2 − σ 2 water ¼ σ 2 − ðk Rrs R rs Þ 2 with the results of comparison for 551 nm shown in Fig. 9 together with k Rrs for both sensors; a reasonable similarity with the smallest σ sag in Hawaii and the largest in the LISCO area is demonstrated. A very small k R rs in the open ocean and an approximately linear increase with GSD in coastal waters are consistent for both sensors.
Examples of OLI processing in the Hawaii, WaveCIS, and LISCO areas are shown in Figs. 10-12, respectively. In addition to the mean R rs ðλÞ (first column), σ spat ðλÞ spectra are shown for the spatial resolution 90 to 5250 m (column 2) and CVs were calculated as a function of GSD and spectrally (columns 3 and 4). OLI observes scenes almost vertically, so the dependence of GSD on the viewing angle is small and was not considered. Three rows in each figure correspond to wind speed ranges W < 3 m∕s, 3 < W < 5 m∕s, and W > 5 m∕s. CVs are presented in two forms: as σ spat ðλÞ∕R rs ðλÞ, which includes all effects, and as k R rs ðλÞ ¼ σ water ðλÞ∕R rs ðλÞ, which is directly related to the water variability.  Fig. 10 Hawaii-near MOBY site, R r s spectra (gray) with their mean values for all GSD (column 1), σ spat ðλÞ spectra for different spatial resolution (column 2), CV ¼ f ðGSDÞ relationship by spectral band and k R rs (column 3), and CV spectra for different spatial resolution (column 4). Data were analyzed using three wind speed (W ) intervals: W < 3 m∕s (row 1), 3 < W < 5 m∕s (row 2), and W > 5 m∕s (row 3).
Very clear waters in the Hawaii area are represented by the station near the MOBY site with spectra shown in Fig. 10 with specific shapes of σ spat ðλÞ, small CVs, and changes in CV values occurring at low GSD below ∼500 m, which then become constant for larger GSD. At the same time, k R rs ðλÞ at this station, as was shown above, is almost equal to zero, meaning that all CVs are related to effects other than water variability.
In the WaveCIS coastal area (Fig. 11), the impact of water variability is much stronger than in the open ocean waters, and k R rs is in a similar range as CVs with some differences depending on the wavelength and wind speed. In the LISCO area (Fig. 12), surface and glint effects are much more pronounced than at WaveCIS (see Fig. 5), and k R rs is usually smaller than CV for all wavelegths. The Moses et al. 22 Fig. 12 Same as in Fig. 10 but for the LIS-station 1. absorption at 440 nm is close to the CV dependence on GSD for the blue band (at least for low and high wind speeds) and is higher than k R rs .

Conclusions
A model was developed to estimate the contribution of the main components that play a role in atmospheric correction, to the variance spectra of remote sensing reflectance σ 2 ðλÞ in different water types. The model was applied to determine the spectral structure of remote sensing reflectance variances due to the spatial variability σ 2 spat ðλÞ and their dependence on GSD based on satellite imagery from SNPP-VIIRS and Landsat 8-OLI sensors. It was shown that, in the open ocean, there was practically no water variability up to GSD of about 1500 m, it remained low for larger GSD, and σ spat ðλÞ was found to be governed mostly by surface and aerosol components. In coastal waters, on average, there were similar variance contributions from surface and aerosol components, whereas water variability played a dominant role and increased proportionally to GSD. The coefficient, k R rs , which characterizes water variability, was typically smaller than CV determined directly as σ spat ðλÞ∕R rs ðλÞ, and therefore water variability cannot be accurately determined from the latter ratio, which is due to the full spatial variability σ spat ðλÞ. Relationships shown in Figs. 8 and 9 allow for estimating these differences. The differences are especially pronounced in the blue bands, where the contributions of surface and aerosol effects are larger.