Spatial Variability of Suspended Sediments in San Francisco Bay, California

: Understanding spatial variability of water quality in estuary systems is important for making monitoring decisions and designing sampling strategies. In San Francisco Bay, the largest estuary system on the west coast of North America, tracking the concentration of suspended materials in water is largely limited to point measurements with the assumption that each point is representative of its surrounding area. Strategies using remote sensing can expand monitoring efforts and provide a more complete view of spatial patterns and variability. In this study, we (1) quantify spatial variability in suspended particulate matter (SPM) concentrations at different spatial scales to contextualize current in-water point sampling and (2) demonstrate the potential of satellite and shipboard remote sensing to supplement current monitoring methods in San Francisco Bay. We collected radiometric data from the bow of a research vessel on three dates in 2019 corresponding to satellite overpasses by Sentinel-2, and used established algorithms to retrieve SPM concentrations. These more spatially comprehensive data identiﬁed features that are not picked up by current point sampling. This prompted us to examine how much variability exists at spatial scales between 20 m and 10 km in San Francisco Bay using 10 m resolution Sentinel-2 imagery. We found 23–80% variability in SPM at the 5 km scale (the scale at which point sampling occurs), demonstrating the risk in assuming limited point sampling is representative of a 5 km area. In addition, current monitoring takes place along a transect within the Bay’s main shipping channel, which we show underestimates the spatial variance of the full bay. Our results suggest that spatial structure and spatial variability in the Bay change seasonally based on freshwater inﬂow to the Bay, tidal state, and wind speed. We recommend monitoring programs take this into account when designing sampling strategies, and that end-users account for the inherent spatial uncertainty associated with the resolution at which data are collected. This analysis also highlights the applicability of remotely sensed data to augment traditional sampling strategies. In sum, this study presents ways to supplement water quality monitoring using remote sensing, and uses satellite imagery to make recommendations for future sampling strategies.


Introduction
San Francisco (SF) Bay is the most extensive estuary system on the west coast of North America, draining approximately 40% of California's land area into the Pacific Ocean [1]. The bay is home to over 500 species of wildlife, and the surrounding area is home to over 7 million people and significant tourism, agriculture, and tech economies. Consequently, monitoring water quality is a vital aspect of Bay-wide management.
SF Bay's connection to a 153,000 km 2 watershed [1], tidal influence from the Pacific ocean, and bathymetry make it a highly dynamic system with continuously evolving water masses. This physical mixing is tied to the distribution of suspended particles in SF Bay [2]. The Delta is the main source of SPM into the Bay [3]. SPM moves with tidal pumping [4], and wind-driven re-suspension over bathymetry drives localized SPM concentrations [5,6]. SPM can also be used as a proxy for a variety of water quality applications: One way to expand SPM monitoring efforts is with remote sensing, which provides broad spatial coverage and the ability to collect a large amount of data regularly. Data collected from ship-based systems (i.e., [16]), airborne sensors (i.e., [9]), and satellite platforms (i.e., [2,8,17], etc.) have been used widely to map SPM. Mapping SPM requires an applicable sensor, as well as a functional retrieval algorithm for converting reflectance data to SPM values. These algorithms are broadly classified as either analytical, semianalytical, or empirical [17,18], where analytical algorithms are physically derived based on optical properties, empirical algorithms are constructed using in situ or in-water data, and semi-analytical algorithms are a combination.
Most remote sensing in SF Bay has utilized satellite remote sensing data and semianalytical or empirical retrieval algorithms. The authors of [2] created an empirical relationship between SPM concentration and reflectance from the Advanced Very High-Resolution Radiometer (AVHRR) to map SPM in the full bay at 1 km resolution. This showed the broad potential of SPM remote sensing to track environmental conditions, specifically flooding and re-suspension events in SF Bay. In recent years, remote sensing technology has advanced and made remote sensing of SF Bay in high-resolution possible. Notably, the authors of [9] used the airborne JPL Portable Remote Imaging SpectroMeter (PRISM) sensor to map turbidity and water quality at high (2 m) resolution, demonstrating the potential of collecting high-resolution data for seeing spatial patterns and variability in dynamic regions of SF Bay. Using an empirical partial least squares regression approach, they quantified and mapped turbidity in Grizzly Bay (Figure 1). While their method requires hyperspectral data, the widely used semi-analytical single-band algorithm by [19] has been used to map SPM and associated pollutants at 30 m resolution with the multispectral sensor on Landsat-8 [8].
Higher-resolution sensors like the European Space Agency's 10 m resolution Mul-tiSpectral Instrument (MSI) on-board the Sentinel-2 satellites provide imagery capable of resolving small-scale features across the entire bay. Sentinel-2 data has been used in other estuary systems with the Case-2 Regional Coast Color processor (i.e., [18])-a semianalytical algorithm that utilizes a machine learning and look-up table approach [20]. This sensor and retrieval strategy are untested in SF Bay, but would provide much improved spatial information.
High-resolution remote sensing data can help expand monitoring by being an additional data source, and by providing more spatial information for making sampling decisions. Understanding the spatial distribution of water quality variables is vital for successful sampling designs [13,21]. Knowing how much variability exists at different spatial scales is an important consideration: for example, coastal systems tend to have most of their variability contained below 250 m, but areas offshore have less small-scale variability [22]. In addition, knowing the amount of spatial variability in different areas of a system is highly important. In SF Bay, longitudinal variability has been examined [13] but bay-wide spatial variability is unconstrained.
With this in mind, our study has two main goals: First, to examine the applicability of remote sensing of SPM in SF Bay and examine how different data types contribute to monitoring in terms of accuracy and spatial coverage. The second goal is to quantify spatial variability of SPM in SF Bay, to contextualize sampling strategies and future monitoring. Specifically, we use remote sensing of SPM to address the question: How much spatial variability exists at different spatial scales in SF Bay, and how much is captured by monitoring at different resolutions? Overall, this project seeks to provide information on (1) how remote sensing can directly contribute to monitoring water quality in SF Bay, and (2) the overall variability that exists at different scales in the Bay.

Study Area
SF Bay can be simplified into three main regions: the North Bay, Central Bay, and South Bay (Figure 1). These regions can be approximately defined by the location of bridges crossing the Bay: North Bay is the region of the Bay north of the Richmond Bridge, Central Bay is the region between the Richmond and Bay Bridges, and South Bay is the region south of the Bay Bridge. These regions each experience different physical processes driving water movement and spatial variability.The bay overall is a shallow wetland system, with most of the bay less than 3 m deep ( Figure 1). The exception to this is a dredged shipping channel which runs through the middle of SF Bay from North to South. This shipping channel is the deepest part of SF Bay (depths greater than 10 m) and is dredged regularly to facilitate the movement of ship traffic. This stark difference in bathymetry and higher flow velocities in the shipping channel [14] make it distinctive from the rest of the bay.
The full SF Bay is connected to two major water bodies: the San Francisco Bay Delta to the northeast, and the Pacific Ocean to the west. The San Francisco Bay Delta feeds into east North Bay and is where the Sacramento and San Joaquin rivers combine. Inflow from the Delta accounts for 90% of the total freshwater inflow to the Bay. North Bay is therefore characterized by a strong salinity gradient, and this inflow of river water is the primary control on bay-wide salinity [1]. The Delta inflow is also the primary source of suspended sediments, nutrients, and freshwater phytoplankton to the estuarine system [23]. The remaining 10% of inflow to the Bay comes from runoff and other point sources.
Flow out of SF Bay is primarily controlled by the amount of freshwater inflow from the Delta [23]. The Bay drains into the Pacific Ocean through the Golden Gate (Figure 1), which exerts strong tidal influence on horizontal mixing within SF Bay moving water back and forth over 10's of km. The combination of inflow from the Delta, tidal mixing from the Pacific, varied bathymetry, and wind-driven mixing creates a highly dynamic system. Each type of physical forcing has different levels of influence on North, Central, and South Bay: for example, the strong salinity gradient set up by Delta flow in North Bay creates a very different environment than South Bay. Water retention time in South Bay is on the order of 3-5 months, while in North Bay water moves through as quickly as 2 weeks [23].

Data Overview
To examine how remote sensing methods can supplement existing in-water monitoring methods, we collected four types of data in collaboration with the USGS Water Quality Monitoring Program. We used in situ SPM grab sample measurements and continuous turbidity flowthrough data taken by the USGS in collaboration with UCSC as part of their standard monitoring procedure; continuous radiometric data taken from the bow of the USGS research vessel, the R/V Peterson; and 10 m satellite imagery from the MultiSpectral Imager (MSI) onboard the European Space Agency's Sentinel-2 satellite. In sum, we used discrete SPM grab samples, continuous flowthrough and shipboard radiometry taken along a transect, and spatially complete satellite imagery covering the full bay.
The ground-based data used in this study were collected on 25 April, 4 June, and 23 October 2019 onboard the R/V Peterson. Each cruise began at approximately 6:00 a.m. local time (PDT) from Redwood City harbor and ended between 3 and 5:00 p.m. in Antioch ( Figure 1). Sentinel-2 imagery was always acquired at 11:30 a.m. PDT. Sentinel-2 overpasses occurred on both the April and June dates. No Sentinel overpasses occurred on 23 October, so instead we used imagery that was collected on 22 October.

USGS In-Water Data
Water quality measurements including SPM concentration (mg/L) were collected at regular point stations ( Figure 1) at 2 m depth, according to USGS standard operating procedures. Quality controlled data from the USGS are available for public use (sfbay.wr. usgs.gov, accessed 21 January 2020).
Continuous underway measurements of turbidity were collected every 5 s at 2 m depth throughout each cruise using a Turner Designs Self-Contained Underwater Fluorescence Apparatus (SCUFA). These underway flowthrough data were validated following the standard operating procedure of the USGS San Francisco Bay Water Quality Monitoring program [12]. Briefly: data were filtered to remove bubble noise in the system. Flowthrough data taken while the ship was sampling at designated stations ( Figure 1) were averaged and compared to corresponding discrete samples using linear regression. If the regression line had an R 2 > 0.7, the relationship was accepted and used to convert all the filtered flowthrough data to SPM units. If the R 2 was less than 0.7, the data was geographically split into groups of at least three points. A new regression was fit for each group and used to convert flowthrough data contained in each group to SPM.

Shipboard Radiometry Data
In addition to the regular USGS data, a HOBI Labs HydroRad-3 (HydroRad) hyperspectral radiometer was mounted at the ship's bow and set to record continuously throughout each cruise. Discrete hyperspectral measurements from an Analytical Spectral Device (ASD) FieldSpec HandHeld 2 Spectroradiometer were also taken to corroborate HydroRad data at locations corresponding to simultaneous water quality measurements.
Prior to being deployed, the data quality and signal to noise ratio of the HydroRad was checked against the ASD in stable conditions. Onboard the R/V Peterson, the HydroRad recorded continuous radiometry measurements of downwelling irradiance (E d ), sky radiance (L sky ), and total radiance from the water (L T ). The HydroRad's three sensors were fixed 2 m above the main deck on the ship bow. Care was taken throughout the cruise to keep sensors pointing 100 • -130 • average solar zenith angle. L T was fixed 40 • down from horizontal, L d was fixed 40 • up from horizontal, and E d was fixed vertically.
HydroRad data collected before station 24 were removed to eliminate data collected in low-light conditions (i.e., solar elevation below 30 • ). The remaining data were processed following [24]. Remote sensing reflectance (R rs ) was calculated from the HydroRad data using the equation: where λ is wavelength, and ρ sky = 0.1 to correct for skylight reflectance off the water surface [24]. Data were inspected for quality control and manually compared to ASD measurements. Flat spectra, spectra with negative values, and oversaturated spectra were removed, and high frequency noise (attributed to glint) was removed from the remaining data using a running mean with a smoothing value of 2. SPM was calculated from HydroRad R rs along each cruise track using the single-band algorithm developed by [19], following the successful use of the algorithm in SF Bay by [8]. This algorithm uses a single channel to determine SPM based on spectral brightness, and then uses an empirically derived set of coefficients to calculate SPM. We chose to synthesize a band centered at 700 nm by averaging all hyperspectral channels within an 8 nm window to mitigate noise effects; 700 nm was deemed appropriate based on the recommendations in [19], preliminary tests, and an observed increase in sensor noise above 700 nm. Retrieved SPM values from data taken at stations were compared to validated flowthrough SPM for each cruise (Figure 2).

Satellite Image Processing
In addition to shipboard radiometry, we tested the applicability of retrieving SPM with Sentinel-2 using existing algorithms. Because our goal was to test the applicability of existing tools to SF Bay in order to identify the best practice for repeated monitoring, no additional calibration was performed on any of the algorithms. To examine SPM outside the shipping channel, we use imagery from Sentinel-2 overpasses on 25 April, 4 June, and 22 October 2019. Concurrent USGS cruises with radiometry occurred on 25 April and 4 June. No ground data was collected on 22 October, but data were available from a cruise on 23 October. A 1 day difference is outside the temporal window for comparison in SF Bay [8]. Consequently, direct comparisons between ground and satellite data are only made for April and June, and the October data are simply used as a qualitative check.
For each date, two images of SF Bay were acquired sequentially at approximately 11:30 am PDT. Sentinel-2 granule 10SEH covers the northern half of SF Bay, and granule 10SEG covers the southern half. No concurrent imagery was available for the Delta region (granule 10SFH) encompassing station 657 for any of our dates because it is imaged on a different cycle. Imagery was downloaded from the ESA open-access data portal (scihub.copernicus.eu) in two different versions; one set of Level 1C top-of-atmosphere (TOA) reflectance, and a second set processed to the ESA's atmospherically corrected Level 2A bottom-of-atmosphere (BOA) reflectance. Images were mosaicked using the ESA's open-source imaging processing software, SNAP, to create full Level 2A image and a full Level 1C image for each date.
We tested three retrieval methods for mapping SPM using the Sentinel-2 imagery. We first tested the single-band ( [19]) SPM algorithm applied to the Level 2A surface reflectance. Per the requirements of the algorithm, we selected band 4 (665 nm) which is the band farthest into the red that has 10 nm spatial resolution.
Due to known issues with the Level 2A atmospheric correction over coastal waters, we also tested the Case 2 Regional CoastColour (CR2CC) processor in SNAP with Level 1C TOA imagery. The CR2CC processor was developed specifically for remote sensing of water, and simultaneously atmospherically corrects imagery and retrieves optically relevant water quality variables (including SPM concentration). C2RCC uses a set of neural networks to match observed conditions with a library of radiative transfer simulations representative of different water types [20]. It has been specifically developed to work with certain satellite sensors, and therefore was not transferable to the shipboard HydroRad dataset.
All atmospherically corrected reflectance products are compared in Figure 3. Surface reflectance from Sentinel-2 was converted from surface reflectance to R rs by dividing by pi to compare with HydroRad R rs . Processing steps for all SPM retrievals are documented in Table 1. Table 1. Processing steps for each retrieved SPM product. Two data sources were used for retrievals: the surface-level HydroRad and Sentinel-2 MSI. Two atmospheric corrections were evaluated for Sentinel-2 ( Figure 3). For one of those atmospheric corrections (C2RCC), two retrieval algorithms were tested. For each product, the absolute error compared to in-water flowthrough data was calculated for both April and June. To evaluate the validity of each atmospherically corrected product, we compared both (the reflectance from Sentinel-2 processed with C2RCC and the reflectance from the standard Level 2A Sentinel-2 product) to reflectance measured by the HydroRad at the surface ( Figure 3). We also compared the satellite SPM retrievals with the HydroRad SPM retrievals and the in-water SPM measurements. In total, we compared five continuous SPM transects ( Figure 4A,C), including: (1) the USGS in-water flowthrough data, (2) the SPM retrieved from the HydroRad using the single-band [19] algorithm, (3) SPM retrieved using the single-band [19] algorithm and the standard Level 2A Sentinel reflectance product, (4) SPM retrieved using the single-band [19] algorithm and the C2RCC Sentinel reflectance product, and (5) SPM retrieved using just C2RCC. For each retrieval strategy, we calculated the absolute error as compared to the flowthrough data to evaluate overall accuracy ( Figure 4B,D).

Spatial Variability Analysis
We examine spatial variability using the SPM distribution maps created from Sentinel-2 imagery using the Ground Sampling Distance (GSD) vs. Coefficient of Variance (CV) method outlined by [21] and modified slightly in this study. This analysis defines the amount of spatial variability at different scales. The method is similar to a semi-variogram or spatial autocorrelation, but we chose this method over others because it is mathematically similar to an inverted Signal to Noise Ratio or uncertainty metric. This allows us to think of the CV as a measure of how representative a single sample is for an area. As larger and larger areas are considered, the variability increases rapidly at small scales first and then more slowly at larger scales, following a logarithmic relationship. , and retrievals three strategies using Sentinel-2 data. SPM from the HydroRad was calculated using the single-band [19] algorithm (solid orange line). SPM was also retrieved using the Sentinel-2 standard L2A reflectance product and the single-band [19] SPM algorithm (light blue line), Sentinel-2 atmospherically corrected with C2RCC and the single-band [19] SPM algorithm (middle blue line), and Sentinel-2 processed entirely with C2RCC for SPM (dark blue line). Panels (B,D) show the overall absolute error associated with each retrieval method compared to the in situ Flowthrough data for April and June, respectively.
The amount of spatial variability is calculated as a function of increasing distance or area (Ground Sampling Distance, GSD). Spatial variability was quantified as the Coefficient of Variance (CV, [21]), defined as where n is the number of segments the dataset is broken into, k is the number of data points in each segment, x i is the within-segment mean, and σ i is the within-segment biascorrected standard deviation. The CV is essentially a measure of how representative a single measurement is for an area defined by GSD.
The rate at which variability is gained with distance can be quantified as dCV/dGSD. We modeled the relationships between CV and GSD for each dataset using a logarithmic function [21]. The transition point where this rate changes (GSD T ) represents the sampling distance at which variability begins to increase more slowly. Moses et al. used two different methods for finding the transition GSD: the first was the find a log-log intersect, and the second was to identify the GSD where slopes reached the 66th percentile. Our data did not display any clear log-log intersect, and the 66th percentile method was arbitrary and did not transfer to our study. Instead, we defined the GSD T as the GSD where dCV/dGSD = 1 × 10 −4 , i.e., where the slope transitioned through 1% CV per 100 m. This gives us a fixed metric to compare across our datasets and represents a clear transition point in dCV/dGSD.
We conducted three distinct spatial variability comparisons. First, we calculated the CV/GSD relationship for the full Bay on each date to examine how variability changed between months. Second, we compared those full Bay relationships with variability in just the shipping channel. This was done by extracting just SPM values along the USGS cruise transect and repeating the variability analysis using just the data within the shipping channel. Third, we examined the variability in different regions of the Bay by repeating the analysis for the North, Central, and South Bay regions.

Environmental Data
We obtained Delta flow data from the California Natural Resources Agency's Dayflow dataset, which is an algorithm used for tracking historical mean daily flows into SF Bay from the Delta (data.cnra.ca.gov, accessed 21 May 2021). Information on tidal phase and state is from NOAA's Tides and Currents Water Level product for SF Bay (tidesandcurrents. noaa.gov, accessed 21 May 2021). Wind data were taken in situ at each station as part of the radiometry data collection protocol. At the time of each satellite overpass, wind was measured at 7.7 m/s, 0.5 m/s, and 1.3 m/s at USGS station 13, which is where the R/V Peterson was sampling at the time.
Tidal conditions on each date were as follows. In April, SF Bay was transitioning from spring to neap tide conditions, and imagery was taken at a peak high tide of 1.6 m. The June image was taken at the height of spring tide, during a falling tide at 0.6 m with the tide going out at time of satellite overpass. In October conditions were switching from spring to neap tide, and it was a falling tide at 0.6 m with the tide going out at the time of the satellite overpass. All imagery was taken at 11:30 a.m. PDT.
Delta flow on each date followed seasonal patterns. The inflow from the Delta was similar for April and June (rainy season), then much lower in October (dry season). Delta inflow in April was 68,464 cfs (cubic feet per second) and in June was 66,005 cfs. The April cruise occurred following a large peak in flow, and the June cruise occurred at the height of a local peak in delta flow. In October, Delta inflow dropped by a factor of four to 15,845 cfs.

Shipboard Radiometry and Flowthrough Comparison
The relationship between retrieved SPM from the HydroRad and in situ SPM from the flowthrough system was correlated (R 2 0.49) on each cruise. This correlation stayed near a 1:1 relationship but differed from cruise to cruise. HydroRad SPM appears to underestimate the flowthrough SPM (slope = 2.15) in April but were more aligned in June and October (slopes of 1.08 and 0.87, respectively) Figure 2.
The full transects of HydroRad SPM, flowthrough SPM, and in situ grab samples are compared in Figure 2. The average distance of grab sample measurements is 8 km, and the average distance between HydroRad measurements was 95.8 m, excluding data gaps and repeats taken when the boat was paused at stations. Notably, a SPM maximum can be identified in all datasets and corresponds with the location and slope of the salinity gradient Figure 2. In April the region of the SPM maximum and salinity transition point are mid-way in the Carquinez Strait near station 9; in June, both features are present just west of station 9; and in October, the transition point moves eastward into the Sacramento River near station 3.
Smaller notable features also occur in the transects. A significant peak occurs near station 9 for all dates, although the exact location of the peak shifts slightly Figure 2. This persistent turbidity feature in the Carquinez straight is not captured by the grab sample stations. The magnitude of this feature changes, from approximately 40 mg/L in April, up to 120 mg/L in June, and down to 15 mg/L in October. Figure 3 shows the comparison of different atmospherically corrected Sentinel-2 products with corresponding surface-level HydroRad spectra. All data were normalized to 740 nm. The C2RCC R rs product was visibly a better match with the HydroRad spectra than the standard Level 2A reflectance product, which consistently over-estimated the overall brightness of the signal.

Satellite Retrievals
Retrieved SPM products from atmospherically corrected Sentinel-2 data are compared in Figure 4 to the HydroRad SPM product and the in-water flowthrough data. For each retrieval strategy, the absolute error as compared to the flowthrough data was calculated ( Figure 4C,D). The standard Level 2A reflectance product with the Nechad algorithm had the highest error. The C2RCC correction had lower error with both the Nechad algorithm and the C2RCC algorithm: in April, the Nechad algorithm outperformed the C2RCC algorithm, but in June the C2RCC algorithm had lower error. Because of the observed difference between the HydroRad and the flowthrough data in April, and the hypothesized stratification confounding results on that date, we elected to use the Sentinel-2 retrieval that had the lowest error on the June date: the C2RCC atmospheric correction and C2RCC SPM retrievals.

Satellite SPM Results
C2RCC was used to generate maps of SPM on 25 April, 4 June, and 22 October ( Figure 5). General statistics for each date are presented in Table 2  Differences in SPM and spatial structure on each date include the widening of the low-SPM area in the shipping channel through San Pablo Bay in the June image compared to the April image. Changes in SPM along the Delta to North Bay gradient through Grizzly Bay occur as the system shifts from high Delta flow conditions (April and June) low Delta flow (October). Closer inspection of the Carquinez straight and Grizzly Bay also elucidate the structure of the Carquinez straight feature observed in the transect data ( Figure 3): a filament of high SPM runs along the northwestern shore of Grizzly Bay and intrudes into the shipping channel cruise track ( Figure 5). This pattern is most pronounced in the June image. The high-SPM filament moves from the shore in April/June to the middle of the channel in October.

Spatial Variability Analysis
To quantify how much variability exists at different spatial scales in SF Bay, we calculated the relationship between spatial scales defined by the GSD and variability defined by CV. Results are displayed in Figure 6. For all datasets, the relationship between GSD and CV followed the logarithmic relationship described by [21]. We considered the modeled CV as a measure of how representative a single measurement is of a surrounding area of size GSD. The rate of change, dCV/dGSD, was consistently positive for all datasets. We specifically address two metrics: the Transition GSD (GSD T ); and the CV at a GSD of 5 km (CV 5 km). The GSD T is the inflection point where dCV/dGSD switches from rapidly increasing to slowly increasing and indicates the scale (GSD) at which variability begins to level off. The CV 5km is the CV at the 5 km scale, which is the approximate scale of USGS point sampling in the Bay.

Comparison between Dates
April exhibited the most variability across all scales (GSDs) for all the dates sampled. Specifically, the CV 5km decreased by 0.2 from April to June and October (full area CV 5km of 0.51, 0.32, 0.32). June and October had much more similar variability with CV 5km of 0.32 for both dates ( Figure 6A,C,E). This is also reflected in the GSD T metric for each date: April, June, and October, respectively, had GSD T values of 921 m, 461 m, and 421 m, indicating much more heterogeneous conditions in April than the other two months.
This pattern of heterogeneity was not reflected in the shipping channel-only analysis ( Figure 6A,C,E). Although April still had the highest variability, June and October had different levels of variability within the shipping channel. CV 5km within the shipping channel was 0.2, 0.16, and 0.12 for April, June, and October, respectively.
Comparing the shipping channel to the full area of SF Bay ( Figure 6A,C,E), we found that there is approximately twice the amount of variability in the full area than just the shipping channel transect. This approximate doubling of variability is consistent across different regions and dates but does change with scale. At smaller scales (i.e., 100 m), the ratio of area variability to transect variability is approximately 1.5; but above a 1 km scale, the ratio is approximately 2.

Regional Comparison
For all dates, Central Bay had the highest CV at all GSDs, followed by South Bay, then North Bay ( Figure 6B,D,F). Specifically considering the average across all dates, North, Central, and South Bay had an average CV 5km of 0.3, 0.53, and 0.43, respectively. At a smaller scale of GSD = 100 m, North, Central, and South Bay each had an average CV of 0.1, 0.18, 0.13 compared to an average full bay CV of 12%. The dCV/dGSD for South Bay was lower than other regions of SF Bay (illustrated by slower increase in CV/GSD relationship, Figure 6B,D,F), implying less small-scale variability. Variability conditions in Central and South Bay were also more inconsistent than North Bay: between all dates at GSD = 10 km, the CV of Central Bay ranged 0.32-0.89, South Bay ranged is 32-81%, and North Bay ranged 0.33-0.37. This indicates that North Bay has the most consistent CV/GSD relationship and spatial structure of the three regions.
The GSD T values for each region ( Table 3) also indicate that North Bay has less spatial variability than Central and South. The GSD T represents the scale at which spatial patterns begin to become more homogeneous and indicates dominant physical mixing processes. A larger GSD T corresponds to more spatial variability in the environment, and a lower GSD T corresponds to less spatial variability. The average GSD T for each region was 451 m, 894 m, and 698 m for North, Central, and South Bays, respectively. A similar but pattern is represented in the transect GSD T values of 118 m, 138 m, and 174 m for North, Central, and South Bays. Notably, the GSD T of the shipping channel transects are lower than the GSD T of the full areas. Also, the transect GSD T s of regions do not consistently correspond to the area GSD T s. For example, in June, the relative patterns of the transect and area GSD T s are the same: South Bay has the highest GSD T for both the transect and areas, and North Bay has the lowest for both. However, in October, the relative patterns of the transect and area GSD T s are different: the Central Bay transect GSD T is smaller than North or South Bay, but the Central Bay area GSD T is larger than the other two regions.

Comparison of Datasets
Shipboard radiometry retrievals matched in situ measurements better than any of the retrievals from the Sentinel-2 data, for both April and June (Figure 4). Changes in the relationship between HydroRad retrieved SPM and validated flowthrough SPM occurred on each date, notably the April matchup relationship had a slope of 2.15 rather than 1 (Figure 2). The discrepancy in the April data could potentially be due to stratification occurring between the in situ grab samples at 2 m depth and the surface SPM concentration retrieved by the shipboard HydroRad. Intensified stratification from high Delta inflow and switching tides could have contributed to the inconsistency between surface radiometry data and flowthrough data at depth. The shipboard radiometry has comparable data to a flowthrough system and is more convenient in that it could easily be set up on small boats or commercial boats already moving around SF Bay. This concept could be highly relevant in SF Bay if scientists partner with commercial vessels.
Retrievals with Sentinel-2 data were able to capture spatial patterns of SPM, with varying success (Figure 4). Using C2RCC to both atmospherically correct and retrieve SPM had relatively low error and was easy to implement using SNAP. C2RCC corrected reflectance matched best with shipboard radiometry, whereas the standard Level-2A reflectance product had much higher brightness across channels (Figure 3). Although C2RCC retrievals did seem to consistently overestimate SPM concentration in April and June, overall it tracked in situ flowthrough data well (Figure 4). The authors of [18] also tested a suite of methods for SPM retrieval at the mouth of the Mzymta River, and also observed a high correlation between C2RCC retrievals and in situ measurements when SPM concentrations were above 20-25 mg/L. In our study (and broadly in SF Bay), SPM concentration was often above this threshold, especially in April and June. However, in October, SPM concentration dropped below the threshold. Since we had no directly corresponding ground data for October, we were unable to confirm the validity of the C2RCC method in SF Bay at lower (<20-25 mg/L) concentrations. More data are required to address this concern and should be noted for future studies.
All sources of continuous shipboard and satellite data captured more variability and SPM features than point sampling. Additional features are visible in the shipboard radiometry flowthrough, and satellite datasets-notably including a peak just after station 9, corresponding to a known gravity circulation feature coming out of Grizzly Bay [7]. At a larger scale, the datasets capture a regional maximum in SPM that corresponds to the location of the X2 salinity marker (Figures 2 and 4). This matches what has been described in previous work as the North Bay Estuary Turbidity Maximum (ETM) zone [7] which serves as an important habitat indicator for estuarine species and marks the transition from fresh to saltwater. X2 tracks the location of the ETM, which in turn is associated with a region of salinity preferred by estuarine species like Delta Smelt [7], and marks the approximate region of transition between marine and freshwater phytoplankton [25,26]. The location of the ETM is dictated by flooding, which controls SPM supply [2,5]. Wind re-suspends SPM, and tides create oscillations in ETM location, which may explain the land-and seaward shifts in the ETM observed in our data. Future work connecting the location of the ETM, spatial variability of SF Bay, and environmental conditions could help inform seasonal requirements for monitoring programs.
The general trade-offs between the methods for monitoring SPM presented in this study are clear: as spatial coverage increases, overall accuracy is lost. This is a function of the additional sources of error that confound SPM with remote sensing (i.e., the accuracy of the atmospheric correction and retrieval algorithm used). It is important for stakeholders to consider what level of accuracy is required to answer various questions and address various issues. Satellite imagery provides a broad and resolved view of SPM distribution; in contrast, point measurements provide certainty that is virtually irreplaceable. As remote sensing technology and retrieval algorithms continue to improve, the trade-off between accuracy and spatial coverage is becoming easier to reconcile. A major focus now is handling and utilizing the shear amount of data and information remote sensing can provide.
While we primarily use shipboard radiometry to connect in situ data to satellite data, it should be noted that shipboard radiometry is a valuable data source in its own right. Shipboard radiometry systems like the HydroRad used in our study provide hyperspectral data. We did not explicitly use the sensor's hyperspectral capacity in this study, but calls for increasing the collection of hyperspectral data in coastal ocean systems are increasingly widespread [27,28]. Consistent hyperspectral data sources could expand biological monitoring in SF Bay and improve our understanding of chlorophyll spatial variability and phytoplankton functional types [27,29]. At minimum, collecting shipboard radiometry provides hyperspectral data with no atmospheric influence at 100 m resolution that could be used for calibration and validation of future hyperspectral sensors like the upcoming Phytoplankton, Aerosol, Cloud, Ecosystem (PACE) and proposed Surface Biology Geology (SBG) missions. Although this study largely overlooks the advantages of having a hyperspectral data source, other studies have demonstrated the potential for retrieving multiple water quality variables from hyperspectral data [9,30]. Empirical retrieval algorithms created using hyperspectral data in Estuary systems are highly accurate when created with sufficient data [31], so having consistent hyperspectral data collection with a shipboard system could facilitate improved retrievals and additional variables to SPM. Using hyperspectral data could help retrieve a suite of water quality variables, which would be particularly helpful for ongoing monitoring of algae and pollutants in SF Bay.

Spatial Variability Differences between Dates
Changes in the CV/GSD relationships between dates corresponded to changes in environmental conditions. Seasonal changes in Delta flow, tidal state, and wind speed in combination with different tidal phase are well known to influence SPM variability [2,6,32,33], and influence GSD T and CVs across scales.
The high SPM concentrations in April and June corresponded to the high Delta inflow as compared to October, and to different tidal phases. Delta flow has been directly related to the overall SPM concentration in SF Bay [2,33]. The authors of [2] documented SPM plumes from large freshwater flows from the Delta in their analysis of AVHRR-derived SPM. It therefore makes sense that April and June had higher SPM at times of relatively high Delta outflow (68,464 and 66,005 cfs, respectively), compared to the low SPM in October corresponding to low (15,000 cfs) Delta outflow. Tidal phase is also linked to high and low SPM concentration: spring tides correspond to high SPM, and neap tides correspond to low SPM [7]. 25 April occurred just before a spring tide, 4 June occurred at the height of a strong spring tide, and 23 October was the beginning of a neap tide.
Daily changes in tidal height also affect overall SPM concentration, but also affect SPM distribution [2,33], which may explain the high variability in April as compared to June and October. The April image was taken during a high tide of 1.6 m, while both the June and October images were taken at a tide height of 0.55 m. Changes in tidal height between images account for significant differences in spatial distribution of turbidity and SPM [9], and could explain the consistently higher CVs observed in April (Figure 6), especially since variability was particularly high in Central Bay which is the most affected by ocean tides.
Changes in wind stress have been shown to affect SPM distribution by controlling re-suspension [6], especially in shallow areas of SF Bay [2]. While we did see higher SPM concentrations in shallow areas of SF Bay on the date with high wind speeds (April, Figure 5, Table 2), this is highly speculative. Wind direction is also important [6] and was not recorded, so more data would need to be examined to draw specific conclusions about the effects of wind on the dates we examined.
We only consider three dates from 2019, but conditions in SF Bay can change dramatically year to year: for example, Delta flow throughout 2020 was much more similar to flow in October 2019 than April or June 2019. Future research should examine spatial structure and variability over a longer time period, during more tidal states, Delta flow conditions and more extensive wind data. Coupling high-resolution SPM remote sensing and physical models of SF Bay (i.e., [6]) could elucidate the relationships between SPM variability and environmental conditions. We recommend using a Principal Component Analysis or similar method to identify regions and times of high and low SPM variability, which could provide valuable information for in-water sampling design.

Differences between Regions
Central Bay consistently had the highest variability, followed by South Bay and then North Bay ( Figure 6B,D,F), which we suggest is linked to its proximity to the ocean and therefore the influence of tides. Because we did not normalize our CV/GSD relationships, and modified the calculation of GSD T , we cannot quantitatively compare our results to [21]. However, [21] observed a decrease in variability in areas with less physical mixing (i.e., offshore), which matches what we observed (i.e., North and South Bay have less physical mixing and less variability than Central Bay).
The results from the CV/GSD relationships between regions show that North Bay has consistently less spatial variability than Central and South Bay ( Figure 6B,D,F), represented by a lower GSD T for North Bay compared to the other regions. This is initially counterintuitive because North Bay is often thought of as highly dynamic compared to South Bay: water residence times in North Bay are approximately 2 weeks, compared to South Bay water residence times of approximately 3-5 months [1]. However, this relatively low heterogeneity matches what is seen in the satellite imagery. For example, if we consider April as a case example: the GSD T of North Bay is 621 m compared to 1.7 km and 1.2 km for Central and South Bay, respectively. Looking at the SPM map for that day, we see that the main source of variation in North Bay is the distinctive shipping channel. However, most of the North Bay region falls within the smaller range of SPM = 60-145 mg/L, whereas both Central and South Bay have a larger range of SPM = 30-145 mg/L. Central and South Bay also don't have just one single distinctive feature the way North Bay does-instead, there is more widespread heterogeneity across the full area of each region, driving up the GSD T . The strong influence of Delta flow in North Bay results in dynamic movement, but steady spatial structure. In contrast, Central and South Bay water movement is less controlled by Delta flow and instead primarily correlated to tidal phase [3,34], which sets up less consistent structure and results in higher spatial variability. South Bay CV/GSD has a lower slope than other regions, which is especially apparent in the transect data.

Differences between the Shipping Channel and Full Bay
In the full bay vs. shipping channel analysis, the variability of the entire bay was consistently higher than the variability observed in the shipping channel ( Figure 6A,C,E).
Results show there is twice as much variability when considering the full bay compared to just the shipping channel transect. Our results suggest that data taken just within the shipping channel will underestimate the variability that exists in the full bay. This can be attributed to two things: first, the shipping channel is a distinctive region from the rest of the bay: it is potentially more homogeneous than the rest of SF Bay because it is consistently deeper (Figure 1) and experiences less wind-driven re-suspension of SPM [6]. Second, the data in the shipping channel were extracted along a transect and are therefore one dimensional while the full area dataset is two dimensional. This is important because variability in SF Bay is not isotropic. Previous studies that have looked at spatial variability in the cross-channel direction [2,9,32] confirm that spatial variability in SF Bay is anisotropic, with more heterogeneity in cross-channel datasets than within-channel datasets.
Our CV/GSD analysis assumes isotropic conditions in the Bay and overlooks the longitudinal anisotropy that exists because of the shipping channel's distinctiveness. Anisotropy could be examined more effectively using a different spatial variability analysis, such as a semi-variogram method that considers the North-South anisotropy introduced by the shipping channel. This can be done by including the direction of strongest correlation and the anisotropy ratio. Despite its limitations, the CV/GSD analysis does effectively illustrate the potential issues with monitoring the Bay only along the shipping channel transect: the channel is physically distinctive and sometimes less heterogeneous than the full bay. Integrating satellite image processing with current monitoring methods would not only provide spatially resolved information, but also provides context on cross-channel variability that is not regularly tracked.

Implications for Monitoring
Monitoring at different resolutions in different parts of SF Bay could be informed by this analysis. Like results from other coastal areas [21], high variability at sub-kilometer scales is observed in SF Bay. In a broad analysis of coastal oceans, the authors of [35] cautioned that for a 30-1000 m area, a single in situ point measurement can have errors up to 18%. In this study, we show that point measurements at 30-1000 m could have much higher error (20-40%) depending on the conditions at any given time. The highly dynamic nature of the SF Bay estuary system also means that there is essentially no point at which spatial structure becomes uniform. Therefore, monitoring decisions about sampling resolution and strategy should be made on a case-by-case basis, based on what processes are being examined and when and where the monitoring is occurring. This is certainly true for SF Bay, and is likely broadly applicable to other estuaries and highly dynamic systems.
Several recommendations emerge from this study. First and foremost, we emphasize the lack of variability captured in the shipping channel. Spatial variability in SF Bay is not isotropic, and the spatial structure of the shipping channel is not necessarily representative of the spatial variability in the rest of SF Bay. Considering data taken at dispersed stations throughout the area of SF Bay is important for any research looking to assess water quality changes or changes in SF Bay conditions. Reliance on point measurements or transect data taken solely in the shipping channel increases the likelihood of missing important outflow events or seeing peaks that are part of a larger feature or small filament.
Similarly, monitoring in different areas of SF Bay should consider the scales at which processes occur, and consider how much inherent uncertainty is incorporated with monitoring at different resolutions. For example, using the CV/GSD relationship established for April, Central Bay has a CV of 0.1 at 20 m resolution, whereas the equivalent CV in North Bay is much higher at 190 m resolution. Accordingly, monitoring in Central Bay must occur at higher spatial resolution than North Bay in order to capture the same level of variability. While this highlights the advantages of high-resolution datasets-like that of Sentinel-2-and provides information for in-water monitoring strategies, our study does not inform on how sampling points could be optimally distributed. We recommend future studies build off this work by using satellite imagery to identify unique and distinct features at different times of year, different tidal states and different Delta conditions.

Conclusions
Monitoring stations placed 5 km apart diagnose general water quality and keep track of changes in overall Bay conditions but miss multiple features at smaller spatial scales and do not capture the full variability of SF Bay. Remote sensing data from satellites and from ships can supplement existing monitoring programs by providing increased resolution and broader spatial coverage. We demonstrate this by comparing in situ measurements with SPM retrievals from a ship-based radiometery system and Sentinel-2 imagery collected on three dates in 2019. Retrievals from the ship-based radiometry data were more accurate than retrievals from the satellite data, but using the C2RCC processor for both atmospheric correction and SPM estimation we were able to map spatial patterns with Sentinel-2. Images of SPM spatial distributions were used to understand spatial variability in different areas of SF Bay in April, June, and October-each of which was associated with different physical environmental conditions. We found that Central Bay consistently has the most spatial variability across all scales, while North Bay has the lowest heterogeneity. In addition, we show variability in the main shipping channel of SF Bay is approximately half of the variability represented in the full Bay on any given date. We suggest that the strength and type of physical processes occurring in different regions of SF Bay are what account for different levels of heterogeneity, and caution that variability in SF Bay is not isotropic. Funding: Partial funding was provided by NOAA Award NA14NOS0120148 and NASA Award 80NSSC18K0420.
Data Availability Statement: Data collected explicitly for this study are available at oceandatacenter. ucsc.edu/outgoing/Taylor_RS_data.zip. All other data used in this study are publicly available as listed in the text: Sentinel-2 imagery is available from scihub.copernicus.eu (accessed on 21 April 2020); USGS data are available at sfbay.wr.usgs.gov/water-quality-database (accessed on 21 October 2019); Delta flow data are from data.cnra.ca.gov/dataset/dayflow (accessed on 21 May 2021). Tide data are from tidesandcurrents.noaa.gov (accessed on 21 May 2021).