Field data, which were used in bio-optical modeling, comparisons of modeled and field Chl-a and other parameters, included data from several Chesapeake Bay cruises. A very comprehensive dataset was acquired by the CCNY-NOAA group in August 2013 at 43 stations, which included Chl-a, inherent optical properties (IOPs) and reflectance spectra. Attenuation and absorption of water and CDOM spectra were measured by the ac-s instrument; backscattering at 5 wavelengths was measured by the bb-9 instrument, both included in the WETLABS (Philomath, OR) package. At each station, upwelling radiance L_u (λ,0^-) was measured using a fiber bundle placed just beneath the water surface and connected to a GER spectroradiometer (SpectraVista, NY). The downwelling radiance above the surface L_d (λ,0+) was measured by pointing the same probe bundle onto a Spectralon plate and the downwelling irradiance was determined as E_d (λ,0+) = A·π·L_d (λ,0+), where A = 0.99 is the reflectance factor of the Spectralon plate (Labsphere, NH), constant for the spectrum in the range of wavelengths from 400 to 800 nm. The underwater remote sensing reflectance R_rs - is then calculated as L_u (λ,0^-)/E_d (λ,0^-) sr^-1, which was adjusted for the propagation through the water-air interface to calculate above surface R_rs . Chl-a from the samples that were collected during the field campaign were determined according to NASA protocol for fluorometric Chl-a determination (Ocean Optics Protocols, 2003).

Capturing the timing and location of a bloom is difficult, and often missing in routine monitoring datasets. An opportunistic sampling event occurred on May 18, 2021, during a high biomass (reaching up to 50 million cells/L) bloom of Prorocentrum minimum. R_rs and water samples for Chl-a were collected at 5 stations in the Upper Bay. Chl-a concentrations were measured in the range of 73–161 mg/m³ and coincided in time with VIIRS and OLCI overpasses. Additional R_rs spectra and Chl-a were acquired in the Chesapeake Bay sporadically from 2014–2016, capturing a range of Chl-a of 9–48 mg/m³. In all these measurements, R_rs were determined from below water HyperOCR depth profiles.

There was also a large dataset of NCCOS R_rs measurements but without corresponding Chl-a. All four R_rs datasets are shown below in Figure 1 in the discussion of the bio-optical model. Ranges of many parameters, necessary for the model and absorption spectra were taken from the previous Chesapeake Bay cruises (Magnuson et al., 2004) and NASA bio-Optical Marine Algorithm Dataset (NOMAD) database (Werdell and Bailey, 2005). Finally, Chl-a data from the Chesapeake Bay program (https://www.chesapeakebay.net/) at multiple stations were used for the validation of the satellite and in-situ data, where most of Chl-a fell below 20 mg/m³.

Concurrent water samples to extract Chl-a and hyperspectral R_rs were collected throughout LIS in 2018–2022 in collaboration with the Connecticut Department of Environmental Protection (CTDEEP) (Turner et al., 2022; Sherman et al., 2023). Additional data were collected from small boats. Hyperspectral R_rs were measured using a HR512-I spectroradiometer (SpectraVista, NY).

Remote sensing reflectance (R_rs ) for VIIRS and OLCI satellite sensors were assessed through comparisons with SeaPRISM instrument (CIMEL Electronique, France) data at the Chesapeake Bay and Long Island Sound (LISCO) stations, where SeaPRISM radiometers are deployed on offshore fixed platforms and are part of AERONET-OC network (Zibordi et al., 2009, 2021). Normalized water-leaving radiances, nL_w (λ), following AERONET-OC protocols and incorporating BRDF correction based on open ocean approaches (Zibordi et al., 2009, 2021), were acquired from the AERONET-OC website for the designated sites. These radiances were transformed into remote sensing reflectance at specific wavelengths. The Long Island Sound Coastal Observatory (LISCO) site (Harmel et al., 2011) upgraded its sensor head in August 2021 to match OLCI sensors with bands at 412, 443, 490, 510, 560, 620, 667, 681, and 709 nm, for detailed comparisons with OLCI data.

The ocean color data employed in this analysis were derived from version 3 level 1.5 data, which underwent cloud screening and quality control measures to ensure data accuracy. All satellite-to-in situ matchups were conducted within a ±2-hour window around the satellite overpass time (Zibordi et al., 2009, 2021).

To develop the NN algorithm, datasets, which connect Chl-a, IOPs and R_rs (λ), water reflectance spectra were simulated based on the bio-optical model (Gilerson and Huot, 2017) with R_rs (λ) including the sum of elastic R_rs ^e (λ) component and fluorescence component R_rs ^f (λ); the latter was included because it is a part of the reflectance detected by the broad I1 600-680 nm band. R_rs (λ) spectra were simulated with 1 nm resolution in the range of 400–750 nm. The maximum of the peak of the fluorescence emission was assumed at 685 nm, fluorescence quantum yield was assumed 1%; the spectral shape of fluorescence was modeled as a Gaussian spectral profile centered at 685 nm, having a full width at half maximum (FWHM) of 25 nm (Mobley, 1994; Gower et al., 2004).

Above water elastic R_rs ^e (λ) was calculated following Lee et al. (2002):

where R r s − (λ) is the remote sensing reflectance due to elastic scattering just below the surface, which is calculated as:

where a(λ) (m^-1) and b _b(λ) (m^-1) are the total absorption and backscattering coefficient spectra, respectively. Broadly used empirically derived parameters (Lee et al., 2009) g ₁ = 0.125 and g ₂ = 0.089, which work well for moderate open ocean and coastal waters were replaced with g ₁ = 0.23 and g ₂ = 0.089 equivalent to the relationship based on our previous studies for a broader range of water parameters (Gilerson et al., 2007, 2015).

The total spectral absorption coefficient, a ( λ ) , is modeled as

where the water absorption spectrum a w ( λ ) was obtained from (Pope and Fry, 1997).

In coastal waters, a p h ( 443 ) , a g ( 443 ) and a N A P ( 443 ) typically have some correlation (even often weak) with each other (IOCCG, 2006). Based on the data from the NOMAD Chesapeake Bay field campaigns (Gilerson et al., 2015) and a p h * ( 443 ) spectra in the Upper Chesapeake Bay (Magnuson et al., 2004) the following relationships at 443 nm were used in the model:

According to Equation 5, a p h * ( 443 ) gradually decreases with Chl-a and remains constant after 60 mg/m³. a g ( 443 ) mostly followed a p h ( 443 ) and a N A P ( 443 ) increases with Chl-a, but less fast than Chl-a itself.

Chl-a were randomly distributed between 0.5 and 200 mg/m³. The spectral phytoplankton absorption coefficient was obtained by multiplying the Chl-a by a specific absorption coefficient ( a p h * ( λ ) , m² mg^-1),

The choice of a p h * ( λ ) strongly influences the corresponding remote sensing reflectance and the emission of fluorescence and was modeled as the specific phytoplankton absorption coefficient in the Upper Chesapeake Bay (Magnuson et al., 2004), shown in Figure 1F with a gradual decrease with increasing Chl-a consistent with Equation 5.

To simulate natural variability, a p h * ( 443 ) were multiplied by a random number drawn from a normal distribution ( N ( μ , σ 2 ) ) with a mean μ=1 and a variance σ 2 =0.04: X 1 ~ N ( 1 , 0.04 ) . In a similar manner, a g ( 443 ) and a N A P ( 443 ) in Equations 6 and 7 were multiplied by X 2 ~ N ( 1 , 0.09 ) . The ranges of variability here and below were based primarily on the published values from IOCCG (2006), NOMAD and the authors’ data for the Chesapeake Bay (Gilerson et al., 2015).

The spectral absorption coefficients of both CDOM and NAP were modeled as having an exponentially decreasing shape with wavelength and referenced to 443 nm (Bukata et al., 1995; Stramski et al., 2001):

S g was modeled as a normal distribution 0.017 N ( 1 , 0.02 2 ) and S N A P as 0.010 N ( 1 , 0.01 2 ) . Equation 7 was also used to determine the concentration of NAP, [NAP] (g m^-3):

where a N A P * ( 443 ) (m² g^-1) is the mass-specific absorption coefficient of NAP at 443 nm, which was simulated as a uniformly distributed random number 0.03 ≤ a N A P * ( 443 ) ≤ 0.05 (m² g^-1). The [NAP] was typically in the range of 0–30 g m^-3.

The total scattering coefficient ( b ( λ ) , m^-1) was simulated as a sum of three components:

Scattering by NAP was modeled using a power law function (Stramski et al., 2001; Twardowski et al., 2001) as follows:

where b N A P * ( 550 ) = 0.5 N ( 1 , 0.04 ) (m² g^-1) is the mass-specific scattering of non-algal particles at 550 nm, and γ 2 = 0.8 N ( 1 , 0.0049 ) .

The scattering by phytoplankton was calculated as the difference between their attenuation and absorption coefficients (Voss, 1992; Roesler and Boss, 2003):

The attenuation coefficient itself was modeled as a power law function (Voss, 1992),

where c p h ( 550 ) = 0.3 C h l a 0.57 and γ 1 = 0.8 .

In the simulations, the backscattering coefficient ( b b ( λ ) , m^-1) was modeled as the sum of the contributing components,

where b b w ( λ ) is obtained according to Morel, 1974 and b ˜ b   p h and b ˜ b   N A P are backscattering ratios for phytoplankton and non-algal particles assumed to be independent of the wavelength (Twardowski et al., 2001; Sydor and Arnone, 1997). Typical values were used as b ˜ b   p h ( λ ) = 0.006 and b ˜ b   N A P ( λ ) = 0.02 .

120000 different conditions were simulated using this model with 70% used in generation and 30% in testing and validation.

As was discussed above, several field R_rs datasets were available for analysis together with (or without) some measurements of water parameters. Four R_rs sets are shown in Figure 1 with corresponding water parameters; some of these parameters (shown in grey) were not measured directly but estimated using available algorithms. Examples of simulated R_rs spectra are also shown in this figure. It should be noted that there was a relatively small flexibility in the selection of parameters described above, which produce spectra similar to the ones in the bloom areas with typical high CDOM and corresponding low R_rs in the blue, spectral features in green-red and a very strong peak around 700 nm comparable with the peak in the green.

In the model development, R_rs (λ) spectra were supposed to be similar not only to the field spectra in Figure 1 , but there were also supposed to be consistent with the good performance of blue-green algorithms for Chl-a retrievals. This should be true at least in the waters with low to moderate Chl-a and RE10 NIR/red bands algorithm for a broad range of waters and Chl-a concentrations, which were observed previously for the Chesapeake Bay (Gilerson et al., 2015).

The ranges of water parameters in the Chesapeake Bay are Chl-a = 0.06–165 mg/m³, CDOM absorption at 443 nm a_g (443) = 0.015–2.0 m^-1, absorption of non-algal particles a_NAP (443) = 0.001–3.0 m^-1, scattering at 443 nm b(443) = 0.3–40.3 m^-1 with the lowest value typically in the Lower Bay and the highest in the Upper Bay (Magnuson et al., 2004). For LIS Chl-a = 1–25 mg/m³, a_g (440) = 0.012–0.5 m^-1, a_NAP (440) = 0.02–0.42 m^-1, particulate backscattering at 650 nm b_bp = 0.005–0.06 m^-1 with the lowest value in the eastern part of the Sound and the highest in the western part (Aurin et al., 2010). In the model Chl-a values were randomly distributed between 0.5 and 200 mg/m³, a_g(443) were mostly in the range of 0–3 m^-1 with decreasing quantities till 6.5 m^-1 and a_NAP (443) = 0–2.5 m^-1.

Several metrics were used in the evaluation of Chl-a algorithms performance which includes a coefficient of determination R², root mean square error (RMSE), relative error e = RMSE/mean as well as recently suggested metrics (Seegers et al., 2018) mean absolute error

and bias

It should be noted that in some figures Chl-a values are shown in the logarithmic scale, while RMSE and e were calculated based on the linear scale.

In continuation of the approach used by (El-Habashi et al., 2016), their simple one-hidden layer multilayer perceptron (MLP) structure was first applied to a newly developed synthetic dataset, to produce a minimum benchmark against which to improve with the introduction of the VIIRS imaging I1 band to complement the 486, 551 and 671 nm band inputs as well as with modifications to the neural network itself. Variables a_ph (443), a_g (443), a_d (443) and b_b (443) were kept as outputs. Chl-a was determined also as an independent output parameter. Performance results are visible in Table 1 . The introduction of the imaging I1 band immediately provided a large performance boost on all four output parameters. However, changes in the neural network structure with the introduction of more neurons in the single hidden layer and the introduction of Rectified Linear Units (ReLU) as the activation function produced a negligible change in the network performance. Similarly, the introduction of a second hidden layer also produced a negligible change in the network performance, indicating that the simpler neural network utilized in previous studies is already capable of capturing the relationships between inputs and outputs well.

In original tests, the bio-optical model was slightly different from the one described above (specific phytoplankton absorption consisted of the micro- and picoplankton absorption with a weighting factor from Ciotti and Bricaud, 2006). In the final version, R² coefficients were higher as shown in Table 1 in parentheses. Figure 2 contrasts the performance of the NNs in the 3-band and 4-band versions against the expected values for a_ph (443), a_g (443), a_d (443), and b_b (443) as measured during the CCNY 2013 cruise in the Chesapeake Bay. In all cases, including the VIIRS imaging I1 band noticeably improves the retrieval quality. In these tests, Chl-a were determined from a_ph (443). When Chl-a were used directly as one of the retrieval parameters, R² for Chl-a became 0.984.

If large datasets of Chl-a and R_rs are available for relevant water conditions, the training can be carried out directly to retrieve Chl-a and other water parameters from R_rs spectra (Hieronymi et al., 2017; Pahlevan et al., 2020). While only 70 points of the field data were available for the Chesapeake Bay, the training gave results quite similar to the ones from the bio-optical modeling, however, some additional tuning was still required, and this option was not further explored.

A total of 70 matchups were included in the tests (43 from CCNY 2013, 22 from Ondrusek 2014-16, and 5 from Ondrusek 2021 measurements) for the validation of NN3, NN4 and VIIRS standard OC3V algorithms on the field data collected across the Chesapeake Bay. Results are shown in Figures 4A–C . The performance of chlC and RE10 on the same field dataset is shown in Figures 4D–F . Among the first three algorithms in Figure 4 the NN4 algorithm shows better performance, although it is worse than the performance of chlC, where all points were used in the tuning and RE10, for which 709 nm band is not available on VIIRS. In Figures 4D–F chlC is plotted against field Chl-a and against RE10; RE10 against Chl-a is also shown for the comparison. High correlations exist for all comparisons in a broad range of conditions in the Chesapeake Bay, but these relationships are not always valid for other types of waters. RE10 was also considered as OC3 VIIRS if RE10 <10 mg/m³.

RE10 was used with the expression (24), which matches the original version in Gilerson et al., 2010, but does not produce complex numbers at low Chl-a

Further tests were performed on SNPP VIIRS data 2012–2022 (NOAA MSL12 data processing) compared with in-situ data from the Chesapeake Bay program (https://www.chesapeakebay.net) and there were 2021 measurements at 5 locations. Results are shown in Figure 5 . The stray light flag was on, HIGLINT and MODGLINT flags were suspended since they did not change the algorithm performance significantly. Most of the points are in the Chl-a range below 20 mg/m³. However, all algorithms, including the OC3 algorithm, retrieve high Chl-a values reasonably well; good performance of OC3 is most likely due to the specific combination of the water parameters in bloom areas, which is not typical for coastal waters with high Chl-a. The time window between satellite and in-situ measurements was ±4 hours. Based on our studies in the Chesapeake Bay, stricter time limits would reduce the number of points but would not improve statistics.

Here and in the figures below the solid grey line marks the 1:1 relationship, while the upper and lower dashed lines mark the limit of Y = X*2 and Y = X/2, respectively, where Y are predicted values and X are expected values.

Performance of the RE10 algorithm for OLCI sensors was evaluated with NCCOS, EUMETSAT, MSL12 and NASA atmospheric correction by the comparison with VIIRS Chl-a in bloom areas with a very broad range Chl-a from 2 mg/m³ to over 100 mg/m³. Because the RE10 algorithm does not provide accurate retrievals for low Chl-a and the OC4 algorithm for OLCI was found not to be always reliable in the waters of the Chesapeake Bay, comparisons were carried out using the RE10M algorithm, where RE10 was replaced with OC3V Chl-a for Chl-a < 6 mg/m³. It was found that the most consistent matchups between VIIRS and OLCI retrievals come from EUMETSAT and MSL12 processing. Examples of such matchups for the Upper Bay and Potomac River bloom areas are shown in Figure 6 . NN4 versus RE10M shows better results than other algorithms. For low Chl-a < 6 mg/m³, OC3V and chlC matchups with RE10M are along 1:1 line because OC3V retrievals are used in all these cases. Since VIIRS algorithms matchups vs RE10M in EUMETSAT and MSL12 matchups produce similar results, both processing approaches from EUMETSAT and MSL12 were recommended for the combined OLCI product. It should be noted that, according to Mikelsons et al. (2022), EUMETSAT processing is more sensitive to the sun glint, which was shown in our comparisons.

Based on the whole study, NN4 VIIRS and OLCI RE10 algorithms were recommended for the combined VIIRS-OLCI product. RE10 was used in combination with OC4 (with OC4 if RE10< 10 mg/m³ and OC4 <10 mg/m³ or OC4 < 10 mg/m³ and clear water conditions based on the diffuse attenuation coefficient threshold K_d (490) < 0.25 m^-1). Examples of the imagery from both algorithms are shown in Figure 7 for May 18, 2021, when field measurements were also available at 5 locations with the coordinates shown in Table 2 , together with measured Chl-a at these points and retrieved from OC3, chlC, NN3, and NN4 algorithms from VIIRS and RE10 from OLCI. Note that part of the area on the OLCI image is masked because of clouds. As before, a slight overestimation of Chl-a is seen in both images in very turbid waters in the Upper Bay, Delaware Bay, and some tributaries. Adjustment coefficients for chlC, NN3, and NN4 are also given in Table 2 . Relative spectral response (RSR) functions were not taken into account in the NN algorithms development to simplify tuning of the algorithms based on the bio-optical model only and comparison with field measurements; for the same reason RSR for the I1 band was considered as RSR = 1.0 for the whole range of wavelengths 600–680 nm. The actual RSR for this band is close to RSR ≈ 0.9, which matches the adjustment coefficient for NN4. The NN3 and chlC algorithms provided similar images but with some adjustments of coefficients, which were less stable than those from the NN4 algorithm. Other examples of images from VIIRS and OLCI for bloom conditions on May 21, 2021, and non-bloom conditions on April 4, 2024, are shown in Figure 8 .

The distribution of absorption and backscattering coefficients at 443 nm retrieved from NN4 together with the SPM concentration based on Equation 20 for May 18, 2021, are shown in Figure 9 , providing additional information about water parameters in the Chesapeake Bay and specifically in the bloom areas, which helps to understand bloom conditions in more details. As can be expected, a_d(443), b_b(443) and SPM have similar patterns since they are mostly proportional to the concentrations of non-algal particles, a_ph(443) and a_g(443) are high in the bloom areas.

The NN4 algorithm was developed based on SNPP VIIRS bands, VIIRS on NOAA-20 has several slightly different bands as was shown above, specifically for the NN algorithm there are M3–M5 bands centered at 489, 556, and 667 nm, and I1 band centered at 642 nm and NN4 algorithm required additional tuning. While the effects of spectral differences between VIIRS-SNPP and VIIRS-NOAA-20 at the blue bands are negligible (e.g., within ~0.1% at M2 band), there are large differences at M4 (green) and M5 (red) bands (e.g., ~16% at M4 band for open oceans) (Wang et al., 2020). Over coastal regions, there are important effects of M4 band difference between VIIRS-SNPP and VIIRS-NOAA-20, because R_rs from NOAA-20 (at 556 nm) is usually much closer to the R_rs peak than that from SNPP (at 551 nm). The same NN4 algorithm was used for VIIRS NOAA-20 bands but with the coefficient 0.65. Images for OC3V and NN4 for NOAA-20 VIIRS are shown in Figure 10 and Chl-a are added to Table 2 . Chl-a from OC3V and NN4 in the same scale looks similar to SNPP Chl-a distributions. For OC3V Chl-a values at in-situ measured points match less accurately with a strong overestimation at two points and were not processed at three other points. NN4 for NOAA-20 is less accurate than for SNPP but can be also recommended for the joint product.

The performance of the NN4 algorithm was validated on the field data in Long Island Sound. Field data were acquired during cruises in 2018–2023 and included radiometric measurements and Chl-a (Sherman et al., 2023). Results for different algorithms are shown in Figure 11 . Most of Chl-a values are below 25 mg/m³, the range that was not the main focus of the NN4 algorithm. The best performing algorithm is OC4 followed by OC3 and NN4. However, all these algorithms perform quite well for Chl-a > 2 mg/m³ and much worse below this value. The NN4 algorithm was used with a coefficient of 0.6, while it was 0.9 in the Chesapeake Bay for SNPP VIIRS. The difference in coefficients might be explained by differences in a p h * ( λ ) with shifts in phytoplankton species (including size and Chl-a packaging), between the time periods in LIS and in the Chesapeake Bay. Optical differences in the water may also influence the bio-optical model. More details about this difference should be further studied.

There were few matchups with VIIRS for field data used in Figure 11 . Sherman et al. (2023) had OLCI retrievals corrected with the Polymer atmospheric correction algorithm and a bio-optical model for moderately turbid waters (Steinmetz et al., 2011), which resulted in good agreement with field observation across the Sound. The performance of algorithms was evaluated at the LISCO site for the period of August 2021–May 2022. The SeaPRISM instrument has bands similar to OLCI bands, and there were no direct field Chl-a measurements. Chl-a were estimated by the OC4 algorithm and compared with those from the NN4 and OC3 algorithms with VIIRS bands, with R_rs determined from the SeaPRISM bands using an adjustment based on the relationship between bands from the synthetic dataset. The NN4 and OC3 algorithms perform very consistently in the whole range of Chl-a from 2–25 mg/m³ as shown in Figure 12 . However, there were no in-situ Chl-a data to confirm these retrievals. Images of Chl-a in LIS based on OC3 and NN4 retrievals for VIIRS and OC4 for OLCI are shown in Figure 13 and are very consistent with each other generally confirming the good performance of algorithms in Figure 12 . Both NN4 and OC3 algorithms for VIIRS can be recommended for the joint product with OLCI OC4.