## Materials and Methods
### Composite Processing
As a result, a variety of constraints were applied to ensure stability in the presence of sunglint. Rather than a least square fit, which is sensitive to outliers and insensitive to small changes in dark features, the absolute value (L1 norm) of percentage difference between co-located, non-cloudy pixels of the reference data and scene were minimized. Pixels over land were excluded during the fitting to ensure consistency over water at the expense of accuracy and consistency over land. To avoid overfitting of sunglint, clouds, or surf, the SBAFs are constrained such that they always preserve a reflectance of 1.0 in the input data and cannot have an unrealistic slope. Furthermore, major changes in band ratios were minimized by including a metric that measures changes in band ratios during the normalization, in addition to the misfit to the reference data. Therefore, the resulting SBAFs avoid overfitting sunglint or waves and are constrained to produce physically plausible models in all cases. The MODIS reference dataset used for normalization is a darkest-pixel composite made from MODIS MOD9A1 surface reflectance data [41]. The MODIS data were composited based on the 20 th percentile of brightness in each band for cloud free pixels collected between 1 January 2014 and 31 January 2019. MODIS was chosen over Landsat 8 or Sentinel-2 as a reference dataset for several reasons: (a) the wide swath width of MODIS reduces systematic row/path artifacts compared to Landsat 8 and Sentinel-2; (b) the daily cadence of MODIS allows for fully cloud-free reference composite; and (c) Landsat 8 and Sentinel-2 surface reflectance data produced with LaSRC have significant across-track artifacts over open water [42]. However, MOD9A1 allows negative reflectance as valid data values and often produces negative reflectance over open water due to slight over-corrections. The MODIS reference composite tends to select over-corrected pixels, and as a result, has uniformly negative surface reflectance values in most bands over water. To compensate for this, a constant shift of 0.01 reflectance was applied in each band to the MODIS reference composite such that the lowest possible value represented in MOD09A1 has a reflectance of 0 instead of \(-\)0.01. This preserves relative color information that is lost over all water pixels if reflectance values below zero are clipped to 0. Normalization reduces major differences in scene-to-scene brightness, but small differences between scenes may create scene edges (Figure 3b). Consequently, a seamline removal method was applied to reduce the effect of scene edges on the final composite (Figure 3c) [43]. Traditional blending or feathering methods operate on a spatially varying weighted average between multiple scenes. However, due to the rapidly varying nature of the ocean surface, as well as the \(\sim\)10-m geolocation accuracy of PS imagery, blending different scenes often introduces more artifacts than it removes. Therefore, seamlineswere removed via a surface reconstruction method that never blends or averages pixels from different scenes, but rather, it smoothly shifts values near scene edges such that the values are equal at the boundary between scenes. This is accomplished with a variant of Poisson surface reconstruction [44], where the gradient of the composite is set to 0 along scene boundaries, and the Poisson's equation is solved for values that honor the modified gradient, while preserving the original values at the edges of each composite tile. While this approach avoids blending multiple scenes, it can shift band ratios significantly in some scenarios. Seamline removal is applied after normalization; however, the changes are minimal, and therefore, major changes in band ratios are avoided. When differences between adjacent scenes are large (e.g., due to unslint, clouds, or surf), seamline removal is only applied when the values along a scene edge differ by less than 0.1 in reflectance.