Spatial Variability of In Situ Above-Water Reflectance in Coastal Dynamic Waters: Implications for Satellite Match-Up Analysis

The validation of ocean color satellite retrievals generally relies on analyzing match-ups between in situ measurements and satellite retrievals. These analyses focus on the quality of the satellite data, however, of the same importance is the quality of the in situ data. Here, we present the spatial variability of in situ above-water reflectance (Rrs(0+)) within the spatial resolution of different ocean color satellites—300, 900, 1500, and 3000 m spatial resolutions, mimicking Sentinel 3 OLCI and MODIS-Aqua satellites, and possible 3 × 3 and 5 × 5 windows. Radiometric data was acquired with autonomous radiometric sensors installed in the British Columbia Ferry Services Inc. vessel “Queen of Alberni” from May to September 2019, crossing the optically dynamic waters of the Strait of Georgia, Canada. The dataset followed optimal geometry of acquisition and processing, including corrections for skylight radiance signals, ship superstructure, the non-isotropic distribution of the water-leaving radiances, and quality control. A total of 33,073 spectra at full resolution, corresponding to 10 days, were considered for the analysis presented here. The results showed that, overall, the subpixel variability increased as the spatial resolution of the sensor or the window size increased, mainly in a linear fashion. Specifically, spatial variability of Rrs(0+) was the largest (∼18% and 68% for 900 and 3000 m pixel resolution, respectively) in Near Field Plume Interface waters, followed by in the Ocean Water Interface (∼28% and 35%, respectively), thus indicating spatial heterogeneity of interface waters. Further, we found that the estuarine waters showed higher subpixel Rrs(0+) variability (∼8% and 16% for 900 and 3000 m, respectively) compared with plume and oceanic waters. We showed that the high spatial variability in Rrs(0+) was primarily associated with the spatial dynamics of the optical water constituents, thus limiting the use of these datasets as Fiducial Reference Measurements and for validation of satellite-derived atmospherically corrected reflectance. We suggest that spatial variability of the in situ Rrs(0+) should also be considered in the selection criteria for good match-up data, especially for data acquired in coastal dynamic systems. As a result, it will advocate for the exclusion of interface or transition water pixel grids in order to avoid compromising the statistical result of satellite validation.


INTRODUCTION
The coastal oceans are highly impacted by climate change stressors, including sea-level rise, warming oceans, ocean acidification, and human activities, leading to eutrophication in coastal regions (Haines et al., 2006;He and Silliman, 2019). In recent decades, the importance and requirement of coastal ocean monitoring systems have been well recognized for the understanding of coastal dynamics and management. Ocean color satellite observations have allowed unparalleled synoptic observations of the biogeochemical variability of the ocean, allowing for a better understanding of ocean dynamics and management (e.g., Boss et al., 2004;Boss et al., 2019;Woźniak and Stramski, 2004;Loisel et al., 2006;Platt et al., 2008;Moore et al., 2009;Dierssen and Randolph, 2013;Sathyendranath et al., 2017;Dutkiewicz et al., 2019). However, the effectiveness of the satellite-derived observations, especially for dynamic coastal waters, requires validation of the different retrievals, including atmospherically corrected reflectance values and reflectancebased biogeochemical variables, such as chlorophyll concentration (Chla) and total suspended matter (TSM) (e.g., Moore et al., 1999;Nechad et al., 2010;Blondeau-Patissier et al., 2014;Shanmugam et al., 2018;Nasiha et al., 2019;Balasubramanian et al., 2020;Giannini et al., 2021;Tilstone et al., 2021). Coastal and estuarine waters are optically more complex as a result of the additional influences of river discharge rich in terrestrial suspended particulates, re-suspended sediments, Colored Dissolved Organic Matter (CDOM), which vary independently of the phytoplankton assemblage (IOCCG, 2000;Davis et al., 2007;Mélin et al., 2007;Loos and Costa, 2010;Blondeau-Patissier et al., 2014;Loos et al., 2017;Phillips and Costa, 2017).
Given this optical complexity, validation is the process of determining the spatial and temporal error fields of a given biological or geophysical data product (Bailey and Werdell, 2006). The satellite match-ups are essential for proper vicarious calibration of sensors, atmospheric correction of top of atmosphere measured radiance or reflectance, and, consequently, accurate retrievals of biogeochemical products (Zibordi et al., 2015a;Zibordi et al., 2015b;Tilstone et al., 2021). In general, these analyses focus on the quality of satellite data; however, the quality of the in situ data is equally important (Tilstone et al., 2021). The quality of ocean color products depends on many factors influenced by various artifacts such as cloud shadows, terrestrial aerosols, exceptional plankton blooms, shallow water with bottom reflection (IOCCG, 2019). Similarly, in situ measurements may have inherent uncertainties including spatial variability, temporal mismatch of satellite overpass, the impact of waves (Zaneveld et al., 2001;Zibordi et al., 2004;D'Alimonte et al., 2013), sensor tilt and self-shading (Gordon and Ding, 1992;Zibordi and Ferrari, 1995;Shang et al., 2017), ship perturbations due to a deployment superstructure (Doyle and Zibordi, 2002;Wang and Costa, 2018;Wang and Costa, 2022), and data processing considering bidirectional effects [Bidirectional Reflectance Distribution Function (BRDF)] of radiometric quantities as a function of illumination and viewing geometries, among others (Park and Ruddick, 2005;Lee et al., 2011;Wang and Costa, 2022).
In situ R rs (0+) measurements are estimated from radiometric parameters such as downwelling irradiance (E d ), sky radiance (L i ), and upward radiance (L t ) and followed by wind correction (Mobley, 1999). For R rs (0+) data satellite match-up analysis, E d , L i, and L t measurements from in situ point stations, located approximately at the centre of the satellite image pixel, are used for validation. The lower temporal and spatial resolution of the imagery acquisition may lead to mismatches between in situ and satellite measured data, especially in coastal waters with riverine inputs where spatial heterogeneity may happen at fine spatial scales (Moses et al., 2016). The uncertainty sources linked to in situ measurements are associated with the spatial representativeness of the in situ observation within the spatial resolution of the imagery (Mélin and Franz, 2014;Pahlevan et al., 2016). For coastal waters, the spatial scales of the variability of the in situ measurements and the variability within the pixel resolution (or average window used) of the satellite imagery have to be considered (Mahadevan and Campbell, 2002;Moses et al., 2016;Barnes et al., 2019).
This paper aims at addressing the R rs (0+) spatial variability of coastal waters within the spatial resolution of different ocean color satellites. Here we evaluate an in situ R rs (0+) data set (N = 33073) acquired with an autonomous set of radiometers (Satlantic Inc./Sea-Bird-Solar Tracker, SAS-ST) in the Strait of Georgia (SoG) waters, Canada, at 300, 900, 1500, and 3000 m spatial resolutions, mimicking Sentinel 3 OLCI and MODIS-Aqua satellites, and possible 3 × 3 and 5 × 5 windows. We look at 1) the spatial variability of in situ R rs (0+) within transitional water boundaries and 2) the spatial variability of R rs (0+) of different water masses in SoG. Estimates of spatial variability in transitional boundaries and water masses provide quantitative guidance for defining the criterion for using in situ R rs (0+) for satellite data match-up analysis in ocean color research. Our premise is that highly varying optical constituents at transitional boundaries between river plume, salt wedge estuary, and oceanic waters lead to larger spatial variability of in situ R rs (0+), contributing significantly to the differences observed between satellite and field data. This may have significant consequences on uncertainties of satellite match-up analysis, which is used to validate atmospheric corrected R rs (0+) and vicarious calibration. Quantifying the uncertainty contributions from spatial variability of coastal water heterogeneity will assist in increasing the quality of data used for satellite validation purposes.

Study Area
The Strait of Georgia (SoG) is a constituent part of the Salish Sea, between Vancouver Island and the extreme southwestern mainland coast of British Columbia, Canada ( Figure 1). It is a narrow passage, measuring 222 km long and 28 km wide on average, in the north-eastern Pacific Ocean (Harrison et al., 1983;Masson and Peña, 2009). The average depth is 155 m, whereas the maximum depth is approximately 420 m. The Fraser River is the primary freshwater input to the central basin of the Strait, with a maximum discharge in June exceeding 11,000 m 3 /s (Thomson, 1981), and the peak discharge period being from June to July (Pawlowicz et al., 2017). Snow-melt constitutes about twothirds of the total runoff, with heavy terrestrial input, which leads to density stratification in the basin region (Johannessen et al., 2006;Loos and Costa, 2010;Halverson and Pawlowicz, 2011). The river runoff creates an estuarine water circulation with strong salinity entrainment in the region that produces a net transport of salt from the saline ocean water to the overlying brackish layer. The horizontal salinity distribution varies greatly depending on the season and distance from the mouth of the Fraser River (Halverson and Pawlowicz, 2011). The general outflows to the Pacific happen through the Juan de Fuca Strait in the south and the Johnstone Strait in the north (Figure 1).
The seasonality and spatial dynamics of the SoG waters are largely influenced by the Fraser River discharge and ocean productivity. Generally, the region has a mean primary productivity of about 280 ± 20 g C m −2 yr −1 (Harrison et al., 1983;Sutton et al., 2013;Johannessen et al., 2021). In the frontal areas of the Fraser River plume, the productivity is still higher due to the reloaded nutrients supplied by the river, and high chlorophyll concentrations have been reported in the waters closer to the Fraser River discharge (Loos and Costa, 2010;Phillips and Costa, 2017). Seasonally, these waters show chlorophyll-a concentrations typically ranging from <1.0 mg m −3 in the winter to >15.0 mg m −3 during spring/ summer and fall blooms (Suchy et al., 2019;Esenkulova et al., 2021).
Optically, the SoG waters are also seasonally highly dynamic as a result of the ocean productivity and high loads of fine inorganic particles and dissolved organic matter discharged into the Strait (Johannessen et al., 2006;Loos and Costa, 2010;Phillips and Costa, 2017). The biogeochemistry of the riverine waters with increasing turbidity levels results in wavelength dependent high scattering in the visible wavelengths (Loos and Costa, 2010;Phillips and Costa, 2017) and, consequently, high reflectance in the longer wavelengths from 550 to 700 nm (Komick et al., 2009;Wu et al., 2014;Carswell et al., 2017;Phillips and Costa, 2017;Giannini et al., 2021). The estuarine circulation causes the mixing of the riverine and oceanic waters, resulting in more optically variable waters with spectrally dependent scattering and absorption (Loos and Costa, 2010;Phillips and Costa, 2017). As a result, relative lower reflectance in the longer 550-700 nm wavelengths is observed due to dominant light absorption. The ocean water exhibits a seasonal dominance of light absorption caused by CDOM and phytoplankton, resulting in the lowest overall reflectance signal relative to riverine and estuarine waters (Carswell et al., 2017;Phillips and Costa, 2017;Giannini et al., 2021).
The cruise transect ( Figure 1) includes oceanic water from the Strait of Georgia and estuarine and plume waters resulting from the Fraser River discharge (Halverson and Pawlowicz, 2011;Loos et al., 2017;Pawlowicz et al., 2017;Travers-Smith et al., 2021). The Fraser River discharge into the Strait forms a seaward extent of freshwater plume mass, bounded by a zone of strong flow convergence resulting in a salinity front. Hence, it develops a brackish surface layer (termed as estuarine waters in this study) between the freshwater plume region and oceanic waters of SoG (Stronach 1977;Thomson 1981;McCabe et al., 2008;Halverson and Pawlowicz, 2011). The estuarine water spread consists of a near-field plume interface (NFPI) and an oceanic water interface (OWI). Entrainment begins at NFPI, leading to significant variability in salinity. The mixing occurs until the oceanic water interface, where salinity variability is lesser than NFPI, in general. These two interfaces explain the transitional water boundaries (Poggioli and Horner-Devine, 2018;Ward et al., 2020). These transitional waters are characterized by highly dynamic bio-physical and hydro-morphologic conditions, resulting in high sub-pixel variability of in situ R rs (0+).

Data Acquisition and Processing
Radiometric data used for this study was acquired from shipborne autonomous radiometric sensors aboard the British Columbia Ferry Services Inc. vessel "Queen of Alberni" from May to September 2019. The data collection was part of the research program "Ferry Ocean Color Observation Systems" (FOCOS). The ship sails at a speed of about 10 m/s from Duke Point (Nanaimo) to Tsawwassen (Vancouver), over a transact length of 45 km ( Figure 1). A set of three hyperspectral radiometers (SAS-ST) were installed on the BC ferry on a customized platform designed by Ocean Networks Canada (ONC). The SAS-ST was installed at 14 m above the water surface on the ferry platform ( Figure 2). The hyperspectral radiometers collect data in the wavelength range of 350-798 nm, whereas the data analysis is presented in the optical range of 400-700 nm. The spectral resolution and spectral sampling interval are 1 nm. The abovewater measurements include total upwelling sea surface radiance (L t (λ)), sky radiance (L i (λ)), and downwelling irradiance (E d (λ)). Considering the speed of the ferry and the frequency of data acquisition for the sensors, we estimated individual measurements at a spatial resolution of 13 m. The sensors were factory calibrated by the manufacturers (Seabird/ Satlantic, Inc.) before deployment, and bi-weekly cleaning of the lenses was conducted.
As per the recommendations from Mobley (1999), the L i (λ) sensor was aligned away from the sun to avoid direct sunglint. The inbuilt motor base of SAS-ST enables the sensor unit to keep at a viewing-sun azimuth, φv, of 120 ± 5°following Hooker and Morel (2003). The L t (λ) was kept at a viewing angle of 45° ( Satlantic, 2016). For the autonomous data collection, we fixed the geometry and acquisition angles based on the ferry's sailing speed and direction (the ferry sails East), the ship heading, the ship structure and shadowing, and the time of day .
The radiometric data was measured, stored on an onboard computer, and transferred to the Oceans Networks Canada server during each ferry sailing. Hence, the near real-time radiometric data collected by SAS-ST can be downloaded from ONC's Oceans 2.0 portal (https://data.oceannetworks.ca/DataSearch) with the data collection date and time. The PySciDON (Python Scientific Framework for Development of Ocean Network applications; Satlantic, 2010;Vandenberg et al., 2017) framework was used to process the data, including specific calibration files associated with the raw data and flags, such as meteorological (cloudy/rainy conditions), sensor rotator angle, and Sun azimuth angle.
The remote sensing reflectance Rrs(0+) was calculated according to Mobley (1999): where Lt(λ) is the total sea-surface radiance, which is measured by the sea viewing sensor, Li(λ) is the sky radiance, and Ed(λ) is the downwelling irradiance. Lt(λ) includes water-leaving radiance (Lw) and a fraction of sky radiance. ρ s is the reflectance coefficient corresponding to the proportion of sky radiance to the water leaving radiance. ρ s depends on variable illumination and sea surface roughness conditions (Mobley 1999). Since the sky glint (ρ s Li(λ)) can have a similar magnitude to Lw, the choice of ρ s significantly influences the accuracy of R rs (0+) estimations (Mobley 1999). In the PySciDON framework (Vandenberg et al., 2017), wind speed data is required as an input to calculate ρ s in the R rs estimation. ρ s was defined considering the local wind speed measured format two Canadian government meteorological stations: longitude from −123.80°to −123.63°corresponding to Entrance Island, while −123.63°to −123.25°corresponds to Sand Heads CS meteorological station (https://climate.weather.gc.ca/historical_data/search_historic_ data_e.html) ( Figure 1). The framework also addresses ship superstructure influence on R rs (0+) calculation. The ship's superstructure (ship wall) influences the above-water radiometry by introducing signal to the radiance field measured by the sea viewing sensor (L t ). The adopted methods for ship superstructure signal contribution to R rs (0+) are fully explained in Giannini et al. (2021) and Wang and Costa (2022), following methods by Hooker and Morel (2003). After the calibration, application of flags and correction factors, the data is binned considering either time or latitude/longitude resolutions.

Ferrybox Data
The salinity, turbidity, chlorophyll and CDOM datasets used in this study were measured in situ by the "Ferrybox" system, which was installed and maintained by Ocean Networks Canada (ONC) (Travers-Smith et al., 2021;Owens et al., 2022). Salinity (PSU) was measured with a SeaBird SBE45 thermosalinograph, Chl-a concentration (mg m −3 ) with a WET Labs ECO Triplet fluorometer, and CDOM fluorescence (ppb) and turbidity (NTU) were acquired with a WETLabs ECO Triplet BBFL2 sensor. The Ferrybox system is located below the main deck, 10 m from the bow, and the sampling depth was approximately 2 m (Halverson and Pawlowicz, 2011;Wang and Costa, 2022).
The Ferrybox data, including turbidity, chlorophyll, CDOM, and salinity, were implemented with ONC's comprehensive process-oriented quality assurance (QA) and product-oriented data quality control (QC) models (Owens et al., 2022). After the pre-deployment testing of oceanographic instruments, automated quality testing was carried out, including QA/QC-related checks, in real-time or delayed, performed via automated quality control procedures while the instrument was deployed. QC of turbidity, chlorophyll and CDOM data includes automatic delayed-mode testing and manual review. The QA/QC test model supports spike detection and gradient steepness tests. For the specific data considered in this study, chlorophyll, CDOM, and salinity datasets were processed with QC flags 1, or 2, representing "data passed all QC tests" and "data probably good," respectively. Among different QC flags, turbidity data is listed with flag 8, which stands for interpolated values, whereas interpolated data exclusively uses clean data (all values have QC flag 1).
Surface current velocity data used in this study were also acquired by ONC using CODAR systems installed at location the Westshore Coal Terminal station (VCOL), Georgina Point(VGPT), Iona (VION), Point Atkinson (VATK).

Data Quality Control
After data acquisition and processing, the entire 2019 dataset was evaluated for quality. The total continuous deployment days of SAS-ST installed on the Queen of Alberni ferry was 116, from May to September 2019. The ferry runs from Duke Point to Tsawwassen daily at local time, 7.45 a.m.-9.45 a.m.; 12.45 p.m.-2.45 p.m.; and 5.45 p.m.-7.45 p.m. The evening trip data collection was excluded from further analysis due to the low sun elevation angle (<30°), and the morning and noon trip data were further evaluated based on the E d and L t data. Figure 3 shows the R rs (0+) for the morning and noon trips on June 16, 2019, as an example along with the data acquisition geometry of the SAS-ST sensor. As shown in Figure 3C for the noon trip data on June 16, 2019 for a specific data point, the rotator angle is at +30°, keeping the sensor plane at 120°from the sun azimuth, whereas the ship's heading is East at 303.5°( Figure 3C). This geometry allows for excellent data quality ( Figure 3A). For the morning trip ( Figure 3D), the rotator plane is at −10°from home orientation, an extreme angle to keep the sensor at 120°azimuth from the Sun. Since the sensor plane could not go beyond the pre-defined backward limit (−10°), the derived R rs (0+) data was erroneous, as shown in Figure 3B, likely measuring the ship structure. Beyond the issues related to the acquisition geometry, the morning trips also exhibited issues with low irradiance and radiance signals (Figure 4). The lower downwelling irradiance (<92 Wm −2 nm −1 at 480 nm) FIGURE 2 | The autonomous radiometer SAS-Solar Tracker is installed on the deck of the ferry, Queen of Alberni. Downwelling irradiance (E d ), sky radiance (L i ) and total radiance (L t ) radiometric sensors are shown. θv = 50°is the sensor viewing zenith angle and θs = 30°is the low sun elevation angle (data acquisition occurs only above θs value).
Frontiers in Remote Sensing | www.frontiersin.org June 2022 | Volume 3 | Article 876748 ( Figure 4B) and L t (<18000 Wm −2 nm −1 Sr −1 ) ( Figure 4D) in the morning data collection due to the lower sun elevation angle (<30°) led to erroneous estimations in the derived R rs (0+) compared to the noon trip ( Figures 4A,C). Given the lower irradiance and radiance signal combined with the less optimal acquisition geometry, data acquired on the morning trips was filtered out, with only noon trip data from~12.30 p.m. to 2.30 p.m. being adequate for further analysis. We selected only days with clear sky conditions for the data acquired during the noon trips. The data set was processed in PySciDON (Section 2.2), and the correction factors were defined as follows: 1) ρ s factor to address skylight radiance signal contributions to the above water signal (Mobley, 1999 2) The ship's superstructure signal contribution to R rs (0+) was defined according to Wang and Costa (2018) and Giannini et al. (2021), following methods by Hooker and Morel (2003). The ship superstructure correction factor for the Queen of Alberni ferry was 0.00005 .
After all the required corrections, the selected data was checked for the tilt angles due to the ferry's stability. The standard deviation in the tilt of the ferry was ±4°, indicating high-quality data. The low tilt angles (pitch and roll) are due to the stable platform, influenced by the large size of the ferry and environmental conditions of the SoG. The final dataset was subjected to BRDF correction to minimize the non-isotropic distribution of the water-leaving radiances in optically complex waters. This followed methods by Lee et al. (2011) and Talone et al. (2018) and was implemented with Python code by Wang and Costa (2022).

Spatial Analysis of R rs (0+) at Different Spatial Grid Sizes
After data reduction for quality control and deriving R rs (0+), the spatial variability of R rs (0+) was evaluated considering the subpixel variability within the pixel resolution of two ocean color satellites, Sentinel 3 OLCI and MODIS Aqua, and commonly used pixel boxes for match-up validation analysis (Mahadevan and Campbell, 2002;Werdell et al., 2007;Barnes et al., 2019;Giannini et al., 2021). A pixel box is generally defined to account for the spatial variability of the biogeochemical information and because the satellite sensor navigation may not be accurate to the pixel (Patt, 2002). As part of a matchup analysis, the median or average (with additional statistical rules, including thresholds for the coefficient of variance) of satellite-derived R rs (0+) (or biogeochemical retrieval) within a pixel box is calculated and evaluated against a co-located in situ measurement of R rs (0+) (or in situ biogeochemical concentration) (e.g., Patt, 2002;Mahadevan and Campbell, 2002;Werdell et al., 2007;Moses et al., 2016;Giannini et al., 2021). As recently highlighted by Barnes et al. (2019), the uncertainties associated with the in situ data should be considered in the validation approach. Here, we use in situ R rs (0+) continuously measured every 2 s from a ferry at a speed of approximately 10 m/s to define the R rs (0+) variability within the pixel grid resolution (PGR) of 300 m (pixel resolution of Sentinel 3 OLCI), 900 m (3 × 3 PGR and similar to MODIS Aqua), 1500 m (5 × 5 PGR), and 3000 m (3 × 3 PGR, commonly used for MODIS Aqua). The radiometric SAS-ST data was used at full resolution without binning during data processing. The cruise track of 45 km was segmented into various PGRs based on the distance between geocoordinates. Every PGR consists of a number of spectra with respect to the pixel window size. For instance, the number of spectra yielded into each PGR was on average as follows: PGR 300, 900, 1500, and 3000 m yielded 15, 44, 73, and 140 individual spectra, respectively. For these PGRs, and considering the length of the ferry run, the total number of pixel windows was approximately 1605, 536, 321, and 161, respectively. For the simulation of PGR, we considered a fixed-width grid size, estimated from the distance between the geographical coordinates of the ferry transact using the cosinehaversine formula (Robusto, 1957). Note that the number of pixel windows for the PGR slightly varies due to the occasional inappropriate viewing geometry of the sensor. For example, on 2 September 2019, we found that no data were recorded for the location between 49.216°; −123.81°to 49.216°; −123.808 due to issues related to the rotator angle of the sensor plane. Given these occasional issues, we excluded the pixel windows with less than 75% of the average number of spectra corresponding to the pixel grid.
After defining the PGRs and individual spectrum associated with each pixel window, the average and coefficient of variance of R rs (0+) for each pixel window were calculated. The coefficient of variation within a PGR (CV Δ ) was calculated as follows: where μ is mean of R rs (0+) spectra within the pixel grid, σ = (xi−μ) 2 N is the standard deviation from the mean, x i is the ith data point within the pixel grid with the total number of data points equal to N.
This study specifically analyzed the R rs (0+) variability associated with transitional water boundaries between the Fraser River plume and estuarine waters and estuarine and oceanic waters. We further extended the analysis to define the spatial variability for different water masses: plume, estuarine, and oceanic, to account for its contribution to the uncertainties associated with in situ R rs (0+) measurements.

Classification of Water Types in SoG
We classified the water masses along the cruise transect into three types: oceanic, estuarine, and plume waters, following a regional salinity-based criterion (Halverson and Pawlowicz, 2011;Travers-Smith et al., 2021), as follows: Here S Ref is the mean salinity calculated for every transect in the region influenced mainly by oceanic waters (Travers-Smith et al., 2021), this corresponds to a transect length of 5 km near Vancouver Island. The estimated S thresh for oceanic waters was defined as above 25 psu, and waters with salinity lower than S thresh were classified as under riverine influence. The riverine waters (S thresh < 25 psu) were further classified into estuarine waters with higher entrainment and mixing (15-24.9 psu), and plume waters (<15 psu).

Rrs(0+) Variability Along the Cruise Transect: Transitional Waters
To illustrate the significant R rs (0+) variability of the transitional waters (NFPI and OWI), we show the variability of R rs (0+), associated with salinity and biogeochemical variables for 2 days, 2 September and 22 August 2019, but omitted the analysis for the other 8 days of the dataset since the results were very similar. For these 2 days, Figures 5, 6 show R rs (0+) spatial variability measured on an entire ferry run between Vancouver Island (Duke Point) and the Vancouver mainland (Tsawwassen) at Band 3, Band 5 and Band 6 of Sentinel OLCI. The central wavelength of Sentinel OLCI sensor bands 3, 5, and 6 is 442.5, 510, and 560 nm, respectively. For 2 September 2019 ( Figure 5), in general, as PGR increased, so did the R rs (0+) within pixel variability, expressed as CV Δ , with the larger variability ranging from about 30% to 70% for 300-3000 m, respectively. Specifically, at 300 m PGR, the mean subpixel variability was about 3.58%, 3.29%, 3.52% (Band 3 (442.5 nm), Band 5 (510 nm) and Band 6 (560 nm) of Sentinel 3-OLCI), except for a few pixel grids. The exceptional large variability was observed for higher PGR in similar pixel locations (~49.066; −123.422 and 49.1692, −123.598); the salinity data revealed the NFPI and OWI waters at these locations, respectively ( Figure 5B). For NFPI, at 300 m PGR, the salinity variability was 36%, and for higher PGRs, the salinity variability was estimated as 26% (for 900 m), 20% (for 1500 m), and 14% (for 3000 m) ( Figure 5B). This large variability in salinity values in NFPI waters indicates the presence of different water masses within the PGR, that is, high spatial heterogeneity, resulting in the highest R rs (0+) CV Δ (Bands 5) for each PGR, 33%, 18%, 56%, and 68%, for 300, 900, 1500, and 3000 m, respectively. For other locations along the track with lower CV Δ , the subpixel salinity variability was comparatively lower than 1%. For OWI, the salinity variability is comparatively lower than that at NFPI, but still considerable. For different PGR, salinity variability in OWI pixel was estimated as 4% (for 300 m), 6% (for 900 m), 7% (for 1500 m), and 6% (for 3000 m). At 300 m PGR, the variability in R rs (0+) for OWI PGR was estimated as 20% (for Bands 3, 5, and 6), that is, slightly lower than the CV Δ for the NFPI waters; for 900, 1500, and 3000 m, R rs (0+) varied by about 28%, 18%, and 35%, respectively.
In these transitional waters, the R rs (0+) variability results from the high spatial variability of the optical water constituents, which increases as PGR increases. Chla, turbidity, and CDOM exhibited greater variability than non-transitional waters ( Figure 5A).
Similar transitional waters and associated variability were identified on August22, 2019 ( Figure 6); however, their locations varied due to tide, current, and river discharge conditions. NFPI was identified at the pixel centered at coordinates 49.039°and −123.341°with a salinity variability of about 11%, corresponding to a subpixel variability in R rs (0+) of approximately 27% (Band 3-442.5 nm), 16% (Band 5-510 nm), and 32% (Band 6-560 nm) for 300 m resolution. For other locations along the track with lower R rs (0+) CV Δ , the subpixel salinity variability was generally <1%, whereas the subpixel salinity variability of the NFPI was as high as 18% (for 900 m PGR), 17% (for 1500 m PGR), and 22% (for 3000 m PGR), associated with R rs (0+) variability of 21%, 29%, and 32%, respectively, at Band 5. For these waters, similar to the NFPI for 2 September, the optical constituent concentrations showed large ranges; turbidity ranged from 3.23 to 8.30 NTU, CDOM from 5.02 to 7.71 ppb, Chla from 2.82 to 0.852 µg L −1 . For the OWI waters, centered at coordinates 49.142°and −123.589°, the lower turbidity ( Figure 6A) indicated an oceanic dominated water mass, which generally shows lower spatial R rs (0+) variability (see next section). At 300 m PGR, the salinity variability was lower than 1% and associated with a lower R rs (0+) variability of 4%. Still, for this day, relatively higher R rs (0+) variability was observed for larger-scale resolutions; for 900, 1500, and 3000 m pixel grid resolution, the subpixel variability was 8%, 11%, and 17%, respectively, for Band 5 ( Figure 6D). Non-transitional waters generally showed similar R rs (0+) variability as OWI, which was associated with lower variability of optical water constituents for each water mass. For example, oceanic waters showed mean values of turbidity, CDOM, and Chla of 2.71 NTU, 4.46 ppb, and 1.75 µg L −1 , respectively, whereas non-transitional waters showed values of 2.82 NTU, 4.69 ppb, and 2.07 µg L −1 , respectively. The mean surface current velocity at the NFPI site (49.039°and −123.341°) on August 22, 2019 was 50 cm s −1 , which is higher than current velocities corresponding to other waters of the measurement transact.
In summary, specifically for the Salish Sea but similar to many coastal waters worldwide, the NFPI generally exhibited larger R rs (0+) variability than OWI. At 900 m (commonly used as a 3 × 3 window for validations of OLCI; Zibordi et al., 2018;Mograne et al., 2019;Kyryliuk and Kratzer, 2019;Giannini et al., 2021) and 3000 m GPR (commonly used as a 3 × 3 window for validations of MODIS-Aqua; Bailey and Werdell, 2006;Carswell et al., 2017;Hilborn and Costa, 2018), the variability associated with the in situ R rs (0+) for NFPI can be higher than 30%. This higher R rs (0+) variability is especially true for the longer wavelengths, illustrated here as Band 6 (of Sentinel 3-OLCI-560 nm), which tends to show relatively higher differences between in situ and satellitederived R rs (0+) when compared with longer visible wavelengths (e.g., Barnes et al., 2019;Tilstone et al., 2021;Zibordi et al., 2021;Giannini et al., 2021). We suggest that the spatial variability of the in situ R rs should also be considered in the selection criteria for good match-up data, especially for data acquired in coastal waters. The spatial variability of in situ data within the pixel of the satellite image is not considered in the present protocols because the data is limited to the pixel size of the satellite (e.g., 300 m × 300 m, 1000 m × 1000 m). Within these areas, we show a large variability of R rs , especially in the transitional and estuarine waters.
In this study, the defined threshold of the coefficient of variance in situ R rs (0+) is 15%, which is within the range suggested in the literature. Generally, the CV threshold values are defined for the pixel box considering the variability of satellite products within the pixel grid box size centered on the location of the in situ measurement, the time difference between in situ and satellite overpass, and the fraction of valid retrievals from the satellite pixel grid box. Mélin and Franz (2014) defined a satellitederived R rs CV threshold in the range of 15%-20%, which is further recommended by IOCCG (2019). In principle, CV thresholds should be considered for both satellite pixel data and in situ data for an improved satellite validation process. For instance, Bailey and Werdell (2006) have defined a Lwn(λ) CV threshold of 15% for a 5 × 5 pixel grid box of Sea-Viewing Wide Field-of View Sensor (SeaWiFS) imagery and a CV threshold of 5% for the in situ Lwn(λ). Similarly, Zibordi et al. (2009), defined a Lwn(λ) CV threshold of 20% for a 3 × 3 pixel grid box of satellite imagery data and, using temporal in situ Lwn(λ) data from the AERONET-OC site (Acqua Alta Oceanographic Tower -AAOT), the authors defined a CV threshold between 5% and 8% for Lwn(λ) measurements. The examples above define the CV thresholds for in situ measurement based only on the temporal variability of these measurements, and do not consider the spatial variability.
In general, the spatial variability of geophysical products is effectively considered for spaceborne remote sensors, whereas the CV of the pixel box is evaluated (Bailey and Werdell, 2006;Zibordi et al., 2009;IOCCG, 2019), while in the case of in situ measurements, the spatial variability is not adequately characterized. This is a result of often insufficient spatially continuous R rs in situ data in relation to the spatial resolution or pixel box of the satellite imagery. This study specifically evaluated the uncertainties associated with the spatial variability of above-water in situ R rs in a dynamic coastal region. The adopted CV threshold value of 15% is inline with the values suggested in the literature for uncertainties associated with temporal variability of in situ data (Zibordi et al., 2009) and satellite pixel (or pixel box) data (Bailey and Werdell, 2006). This present study is based on high spatial resolution data, which allowed us to analyze the in situ spatial variabilities of reflectance measurements. Hence, this study recommendation will allow for enhanced quality of satellite validation procedures. We recommend the consideration of the spatial variability of in situ R rs (0+) measurements and excluding the measurements where the spatial variability, represented by the CV, is greater than 15% within the scale distance of the pixel grid considered for the satellite validation match up analysis.
The inclusion of in situ R rs (0+) for these interface waters may lead to high uncertainties for satellite validation statistics. For example, Giannini et al. (2021) highlighted issues related to the quality of R rs (0+) match-ups (in situ vs. Sentinel-3 OLCI) due to the movement of water masses in this dynamic coastal system. Therefore, for dynamic coastal regions, especially under the influence of river discharge, it is recommended to consider, if available, the spatial variability of salinity as a guideline to define the acceptance of in situ R rs (0+) as a match-up.
Given the difficulties in acquiring continuous R rs in situ data, we also recommend using the salinity threshold of >5% to identify the water masses with potentially significant variability in R rs (0+), which are generally observed in transitional or interface water regions. As expected, we observed that the greater salinity variability in the transitional or interface waters is proportionate to the greater variability in reflectance of the same water mass. Hence, the salinity threshold can be used as a proxy to identify the recommended CV threshold value for R rs (0+) in constructing the matchup data set.
Salinity data is readily available for our sampling strategy because it is collected continuously with the Ferrybox system (Travers-Smith et al., 2021). This is likely the situation with many research cruise ships (e.g., Koponen et al., 2007;Slade et al., 2010;Westberry et al., 2010;Tilstone et al., 2021). However, spatially continuous salinity data may not always be available with other forms of sampling strategy, for instance, R rs (0+) data collected from fixed platforms such as AERONET-OC . We recommend that pixel grids with salinity variability higher than 5% should be considered as a transitional water interface pixel (NFPI or OWI) for our sampling conditions in the Salish Sea. As such, in situ R rs (0+) for these regions will be excluded for satellite match-up analysis (for all considered PGR 300, 900, 1500, and 3000 m) to minimize uncertainties associated with validation of satellite-derived atmospherically corrected R rs (0+).

R rs (0+) Variability for the Different Water Masses: Oceanic, Estuarine, and Plume Waters
After defining the R rs (0+) variability for the transitional waters, NFPI and OWI, we evaluated the spectral homogeneity of the different water masses. For each water type, we calculated the average R rs (0+) variability [denotes as CV Δ hereafter] corresponding to 10 days from May to September (a total of 33073 spectra at full resolution, i.e., without spatial binning) for four distinct pixel grid resolutions of 300, 900, 1500, and 3000 m for Band 3 (442.5 nm), Band 5 (510 nm), and Band 6 (560 nm) of Sentinel 3 OLCI (Figure 7). In general, the scale of the spatial variability differed for each water type, thus showing the distinct nature of these three water masses. Regardless of the pixel grid resolution, the CV Δ values were generally higher for estuarine waters than for oceanic and plume waters. For example, CV Δ was 4%, 9%, and 7% for plume, estuarine, and oceanic, respectively, at 1500 m PGR for Band 5. Similarly, at 3000 m PGR (3 × 3 of MODIS Aqua), CV Δ for estuarine waters was 15%, whereas oceanic and plume waters showed lower variability at Band 5, 9% and 6%, respectively.
The relationship between the different pixel ground resolutions and the R rs (0+) variability within each PGR, defined by CV Δ for each water type was assessed based on the rate of change, d CV Δ /dPGR (defined as Slope value-SL). SL was consistently positive for the three water types, and the relationships were generally linear ( Figure 8). The SoG estuarine waters showed the highest positive slope value (SL = 0.0036), whereas oceanic (SL = 0.0022 for Band 5 and 6, and 0.0034 for Band 3) and plume waters (SL = 0.0009) showed a relatively lower rate of change, thus indicating that for these waters, the PGR approximated the scale of the optical variability of the water. Specifically, for these waters, at 300 m grid resolution, the CV Δ was the lowest (<5% for the three representative bands), but progressively increased with higher PGR, especially for the oceanic waters ( Figure 8A). Estuarine waters showed slightly higher variability, with CV Δ about 2.5 times higher between 300 m (~6%) and 3000 m (~16%) grid resolution. These findings are similar to Moses et al. (2016), who showed a quasi-linear relationship (moderate change in CV Δ ) between pixel spatial resolution (>200 m) and R rs (0+) variability when comparing imagery acquired by different sensors at different spatial resolutions for different types of coastal waters. However, according to the authors, depending on the local spatial scale of the bio-optical and optical constituents, the more pronounced changes in CV Δ happened for PGR lower than 100 m (coastal waters) and 200 m (offshore waters); still, the largest CV Δ were reported for pixel resolution higher than 400 m. Here, we did not investigate the R rs (0+) variability for PGR lower than 300 m since no current ocean color satellite presents finer resolution.  The CV Δ depends on the spatial scale of the ocean dynamics and the optical constituents' concentration in the different waters (e.g., Moses et al., 2016;Taylor and Kudela, 2021). Here, we show the within-transect R rs (0+) variability (CV Δ ) along with the corresponding variability in turbidity, CDOM, Chl-a, and salinity acquired concurrently on 2 September 2019, as an example of our dataset ( Figure 10). For the same day, remote sensing reflectance spectra are shown in Figure 9, for oceanic, estuarine, and plume water types (full resolution data with no spatial binning). On this day, CV Δ in turbidity was the lowest for plume waters (13%) compared with estuarine waters (19%) but similar to oceanic waters (13%); CDOM variability showed a similar pattern with the lowest variability for plume waters (6%), followed by higher variability in oceanic water (9%), and the highest for estuarine (11%); Chl-a variability observed for the plume, oceanic, and estuarine waters was 4%, 2%, and 5%, respectively. Salinity variability observed for the plume, oceanic, and estuarine waters was 6%, 1%, and 13%, respectively, emphasizing the relatively larger variability of estuarine waters. The combined influence of optical water constituents resulted in the highest R rs (0+) CV Δ of 36% for estuarine water mass, whereas for oceanic and plume waters, it was 8% and 8%, respectively ( Figure 10).
The varying ranges of optical constituents were a result of the seasonal dynamics of this region, where during the spring freshet, the Fraser River plume is rich in particulate inorganic and dissolved organic matter, and the dominant phytoplankton bloom happens in the spring and summer seasons (Johannessen et al., 2006;Loos and Costa, 2010;Allen and Wolfe, 2013;Phillips and Costa, 2017;Suchy et al., 2019). It is noticeable that the plume waters were well mixed, demonstrated by the lower variability of the optical water constituents and salinity ( Figure 10). Further, specifically for the plume waters, the high turbidity, which is generally associated with high inorganic particulate scattering (Loos and Costa, 2010;Phillips and Costa, 2017;Giannini et al., 2021), resulted in the highest and least variable R rs (0+). Shelf circulation, the dilution rate of river-borne materials, and transport processes within the plume lead to stratified-shear mixing (Horner-Devine et al., 2015). Horizontal advection, transport of buoyancy, and momentum of Fraser River waters contribute to the well-mixed plume waters. Previous studies reported highly intense turbulence mixing in the near field region of the plume of the Fraser River (Moum et al., 1995;MacDonald and Geyer, 2004). The turbulent kinetic energy dissipation rate is observed to be as high as 10 -4 -10 -3 m 2 s −3 (MacDonald et al., 2007;McCabe et al., 2008;Kilcher et al., 2012;Horner-Devine et al., 2015). Wind and wave forcing also have a vital role in active mixing in the entire plume region (Houghton et al., 2009). Houghton et al. (2009) have shown that the salt fluxes vary from approximately 5 × 10 -5 kg m s −1 during low FIGURE 9 | R rs (0+) spectra measured on 2 September 2019 for (A) oceanic waters (Min; Max; SD; Mean = 0.0019; 0.0029; 0.00017; 0.0021 Sr −1 (B) plume waters (0.0131; 0.0186; 0.0013; 0.0165 Sr −1 , respectively) and (C) estuarine waters (0.0021; 0.0199; 0.0021; 0.0058 Sr −1 , respectively). The descriptive statistics shown are for Sentinel 3 OLCI Band 5 or 510 nm.
FIGURE 10 | Variability (CV) analysis of R rs (0+) at Band 3 (442.5 nm), Band 5 (510 nm) and Band 6 (560 nm), Turbidity, Chlorophyll, CDOM, and salinity parameters for oceanic, estuarine, and plume waters on 2 September 2019 SoG waters. The analysis shown at 300 m high resolution PGR. wind conditions to 3 × 10 -4 kg m s −1 for high winds (12 m s −1 ). For our observation period (May to September, 2019), high wind conditions (4-13 m s −1 ) dominated, causing well-mixing of the plume of water. These plume waters form a strong vertical stratification at around 5-7 m depth, depending on the tide and season (Halverson and Pawlowicz, 2011), resulting in a highly scattered and homogenous surface layer (Loos and Costa, 2010). Hence, the homogenous surface layer of the Fraser River plume waters caused lower variability in turbidity or total suspended matter, which consequently resulted in lower R rs (0+) CV Δ . The oceanic waters in the transect were less influenced by the river plume and resulted in lower R rs (0+) CV Δ .
The estuarine waters showed the highest variability of the optical constituents ( Figure 10) due to salinity intrusion and, consequently, the largest R rs (0+) CV Δ . Since estuarine dynamics and circulation contribute to the transport and dilution of fresh plume waters to saline oceanic waters (Thomson, 1981;Kostaschuk and Atwood, 1990;Hickey et al., 1998), the optical constituents result from both distinct water masses, with absorption and scattering playing a quasi-equal role in light attenuation (Loos and Costa, 2010;Phillips and Costa, 2017). This biogeochemical and, consequently, optical spatial complexity contributed to the largest R rs (0+) CV Δ at any PGR, indicating that the spatial scale of dynamics is shorter than the PGR.
The spatial variability of the in situ R rs (0+) analysis presented here has implications for the satellite data match-up analysis commonly used to validate atmospherically corrected R rs (0+). As per the results obtained, irrespective of the spatial resolution or window box used in the validation process, in situ measurements acquired in estuarine waters may produce high uncertainties for satellite validation statistics.

CONCLUSION
This research provides a detailed and straightforward analysis of the spatial variability of in situ R rs (0+) based on 33,073 measurements from the Salish Sea, west coast of Canada, taken aboard a BC ferry with autonomous radiometers (SAS-ST). The present study brings awareness to the spatial variability of R rs in coastal waters for consideration when defining the quality of in situ data for satellite validation. The comparison of satellite-derived R rs (0+) with the corresponding in situ data is not the objective of this paper. However, our near-future work will involve the detailed satellite validation match-up analysis aspects of this study. The Salish Sea is a highly optically dynamic water system with a strong influence from riverine and oceanic waters, which is common in many coastal regions of the world. This analysis provides a recommendation for accessing the quality of in situ R rs (0+) match-up data based on the spatial variability as complementary to standard recommendations, including the temporal match-up intervals, sensors specifications and calibration, platform perturbation noise, parameterization to address surface Fresnel reflection of the sky radiance and corrections for the bidirectional effect, sky conditions, and data processing methods (Mobley, 1999;Hooker and Morel, 2003;Zibordi et al., 2012;Talone et al., 2018;Vabson et al., 2019a, Vabson et al., 2019bRuddick et al., 2019;Alikas et al., 2020;Tilstone et al., 2021;Wang and Costa, 2022). For specific water masses, the spatial variability of in situ R rs (0+) within the pixel resolution of ocean color satellites or commonly used pixel window sizes may significantly contribute to the uncertainty budget in a match-up analysis. For the simulation of PGR, we considered a fixed-width grid size of 300, 900, 1500, and 3000 m, whereas the distance between the geographical coordinates of the in situ measurement transact was estimated using the cosine-haversine formula (Robusto, 1957). We consider that including a 300 m PGR is sufficient to identify and define the oceanographic features of the region, even with a fixed width window. However, a moving window analysis could contribute to our future studies of very high-resolution data, such as the validation of MSI Sentinel-2, Landsat-8, etc.
Based on a large dataset covering a variety of optical water masses, we found that the spatial variability of R rs (0+) measurements is the largest (~18% and 68% for 900 and 3000 m PGR, respectively) in NFPI waters, followed by in the OWI (~28% and 35%, respectively), thus indicating the spatial heterogeneity of interface waters. Further, the R rs (0+) variability was evaluated for specific water masses beyond the transitional water-plume, estuarine, and oceanic waters-characterized based on a salinity gradient criteria. Overall, the subpixel variability increased as the spatial resolution of the sensor or the window size increased, mainly in a linear fashion. Similar to the findings for the NFPI waters, the estuarine waters showed higher subpixel R rs (0+) variability (~8% and 16% for 900 and 3000 m PGR, respectively) among the three water masses, again indicating its inhomogeneous spatial nature. The high variability in R rs (0+) was primarily associated with the spatial dynamic of the optical water constituents, thus limiting the use of these datasets as Fiducial Reference Measurements (FRMs) and for validation of satellite-derived atmospherically corrected R rs (0+). We suggest considering the spatial variability of in situ R rs (0+) measurement, which represents the small scale environmental changes in coastal waters. Also, we recommend using, the salinity threshold as a proxy to identify the recommended values of R rs (0+) CV Δ in the construction of the match-up data set. Hence, the transitional water boundary pixels, especially those representing NFPI and OWI, and estuarine water masses will be excluded for satellite match-up analysis to avoid compromising the statistical result of satellite validation.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://data. oceannetworks.ca/DataSearch.

AUTHOR CONTRIBUTIONS
HJN was responsible for the data processing, analysis, and manuscript writing. MC and FG were responsible for the project conceptualization, results, discussions, and significant reviews in the manuscript, and ZW was responsible for the development of methods for acquiring and processing the in situ R rs (0+) data.

FUNDING
HN was supported by the Canadian Space Agency (FAST 18FAVICB09) program. The project also had funding from MC from NSERC NCE MEOPAR-Marine Environmental Observation, Prediction, and Response Network, Canadian Foundation for Innovations (CFI), and NSERC Discovery Grant, Canada.