Coastal and inland waters are the most affected by changes in environmental conditions, industrial development, or human intervention, leading to an increase in dissolved and suspended particles like phytoplankton, sediments, and colored dissolved organic matter in the water column (Stumpf et al., 2012; Brown et al., 2015; Pick, 2016; Binding et al., 2021; Jane et al., 2021). Gradually, the interaction of these dissolved and suspended particles with the aquatic environment leads to major changes in biogeochemical parameters (BPs) such as chlorophyll-a (Chla), total suspended solids (TSS), and colored dissolved organic matter (cdom), which subsequently affects the water quality as well as visibility within the water column (Babin et al., 2003). Managing the water quality of inland and coastal waters poses a challenge due to anthropogenic activities, eutrophication, and climate change. Monitoring the water quality for these aquatic ecosystems is essential for sustainable water resource management, providing the foundation on which it is based (Bartram and Ballance, 1996). However, tracking long-term changes in BPs through ship-based measurements poses significant challenges and expenses.

Over the past 25 years, freely available multispectral radiometric images from a series of multispectral ocean color satellite sensors have compiled a vast global observational database of remote sensing reflectance (McClain, 2009). These sensors, deployed on various platforms, include the Environmental satellite (Envisat) Medium Resolution Imaging Spectrometer (MERIS) spanning from 2003 to 2011, the Moderate Resolution Imaging Spectroradiometer (MODIS) aboard Terra (since 1999) and Aqua (since 2002), as well as the Visible Infrared Imaging Radiometer Suite (VIIRS) on the Suomi National Polar-orbiting Partnership (Suomi NPP, since 2011). These heritage sensors, some of which are still operational at the time of writing, have unequivocally broadened our awareness and comprehension of ocean and coastal ecosystems (Cao et al., 2022).

Central to aquatic remote sensing is the concept of apparent and inherent optical properties (AOPs and IOPs, respectively). IOPs include the spectral absorption and scattering properties of water and its constituents (Gordon et al., 1975). The spectral shape of the total absorption coefficient depends on both water and the BPs. The absorbing IOPs can be further divided as a function of these constituents (Werdell et al., 2018). The absorption and scattering coefficients depend solely on the medium or water column properties and are not influenced by the geometry of the light field (Werdell et al., 2018). Moreover, the absorption by BPs is typically close to zero at near-infrared (NIR) wavelengths, making NIR wavelengths indirectly useful for removing atmospheric interference in satellite-captured images (Bailey et al., 2010), unless in very turbid waters. The optically relevant BPs and IOPs can be estimated from the measured multispectral remote sensing reflectance ( R r s λ ), an AOP obtained after the radiometric, geometric, and atmospheric correction of top-of-atmosphere (TOA) satellite measured radiance (Gordon and Wang, 1994). This quantity, R r s (= L w / E d ), is defined as the ratio of water-leaving radiance ( L w ) and downwelling irradiance ( E d ) evaluated just above the water (Mobley, 1999). The spectral shape and magnitude of the R r s are highly sensitive to the processes of absorption and scattering (IOPs) by optically relevant BPs in the water while being relatively insensitive to changes in the ambient light field (IOCCG, 2008; IOCCG, 2009; IOCCG, 2014), making R r s suitable for retrieving the IOPs.

Numerous remote sensing algorithms exist for quantifying BPs and IOPs from R r s captured using ocean color satellite sensors over oceanic and coastal waters. Empirical algorithms primarily rely on single-band linear relationships or multi-band R r s ratios for estimating BPs, such as Chla (O’Reilly et al., 1998; O’Reilly and Werdell, 2019) and TSS (Nechad et al., 2010; Ondrusek et al., 2012; Novoa et al., 2017). The spectra of absorbing IOPs are estimated using exponential functions with the required input parameters, such as absorption coefficients and corresponding slopes of phytoplankton - a p h (Sosik and Mitchell, 1995; Babin et al., 2003; Bricaud et al., 2004; Devred et al., 2022; Wei et al., 2023), non-algal particles - a n a p (Bowers and Binding, 2006), and colored dissolved organic matter - a c d o m (Binding et al., 2008; Shanmugam, 2011; Mannino et al., 2014) at their reference wavelengths. The slopes are obtained from multi-band R r s ratios or functions of BPs. However, these algorithms require substantial datasets to calibrate the underlying models and are primarily applicable to the specific water type for which the model is calibrated. Outside of empirical algorithms, IOPs are estimated using inversion algorithms grounded in statistical modeling, such as Garver Siegel Maritorena – GSM model (Maritorena et al., 2002), Generalised IOP - GIOP model (Werdell et al., 2013) and Linear Matrix Inversion - LMI algorithm (Hoge and Lyon, 1996), utilizing semi-analytical methodologies and prior knowledge of BPs (Wang et al., 2017). Semi-analytical algorithms primarily utilize R r s to deduce IOPs and subsequently estimate BPs by integrating empirical parameters and bio-optical models. Furthermore, the quasi-analytical algorithm - QAA, widely employed over the past decade, represents another prevalent model for IOP inversion (Lee et al., 2002; McKinna et al., 2015; Montes-Hugo and Xie, 2015).

Advanced machine learning models show promise for non-unique inverse problems, such as retrieving BPs and IOPs from R _rs . One such approach is using Bayesian neural networks (BNNs) (Werther et al., 2022). BNNs can retrieve water quality parameters from satellite imagery and have been validated using satellite data. Another proven machine learning method entails using mixture density networks (MDNs) and has been leveraged for the estimation of both BPs and IOPs (O’Shea et al., 2023). MDNs address non-unique inverse problems, like inferring BPs from R r s , by modeling the output as a probability distribution over possible output values. Therefore, the distribution is modeled using a mixture of Gaussians (Bishop, 1994; 1995). Instead of providing the average value of the expected output distribution, like typical multilayer perceptrons (MLPs), MDN provides the full output distribution, enabling the user to study the distribution and choose an appropriate output value, enhancing the overall estimation when the predicted output distribution is multimodal or asymmetric. Further, since the MDNs and BNNs predict output distributions, these distributions can also be used to model the confidence in a specific prediction (Saranathan et al., 2023; Werther et al., 2022) enabling users to disregard predictions on data outside of the training dataset’s range. Specialized MDNs have been developed to estimate specific water quality indicators, BPs and IOPs, such as Secchi disk depth (Maciel et al., 2023), Chla (Pahlevan et al., 2020; 2021b; Smith et al., 2021), TSS (Balasubramanian et al., 2020), phycocyanin (O’Shea et al., 2021; Fickas et al., 2023), and a p h (Pahlevan et al., 2021b) from both multispectral and hyperspectral satellite data, applicable to inland and coastal waters. Furthermore, this architecture has been effectively utilized to concurrently estimate two BPs and an IOP at a single wavelength [Chla, TSS, and a _cdom (440)] (Pahlevan et al., 2022) and tested with multispectral satellite sensor datasets, such as from the Multispectral Instrument (MSI) and Ocean and Land Colour Instrument (OLCI) onboard Sentinel-2, and Sentinel-3, respectively. Recently, the model has been effectively applied for hyperspectral remote sensing by estimating BPs (Chla, TSS, PC) and hyperspectral IOPs ( a p h , a n a p , and a c d o m ) and tested with satellite data from the Hyperspectral Imager for the Coastal Ocean (HICO) and PRecursore IperSpettrale della Missione Applicativa (PRISMA) missions (O’Shea et al., 2023; Lima et al., 2023), indicating its potential application to the recently launched Ocean Color Instrument (OCI) of the Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) mission. In general, due to their demonstrated effectiveness in the non-unique inverse estimation of BPs and IOPs, along with their inherent capacity to gauge uncertainties, MDNs are highly suitable for concurrently retrieving various BPs and IOPs for scientific and water quality management applications.

Daily satellite visits offer a wealth of information crucial for retrospective analysis, particularly in the realm of multispectral ocean color data, which can be effectively used for studying climate change and biogeochemical modeling (Behrenfeld et al., 2006; Platt and Sathyendranath, 2008; Stramski, 2008; Hu et al., 2010; Wang et al., 2011). MODIS has been providing long-term ocean color data since its mission began (1999-), shortly followed by MERIS sensor (2002-). After MERIS’s operational lifeended (2012), the ongoing VIIRS sensor (2011-) continued to provide ocean colordata together with MODIS (Lai et al., 2023). Accurately assessing and quantifying uncertainties in satellite ocean color products is crucial, making thorough validation a challenging task for ocean color missions (Hooker and McClain, 2000). Since these are long-term satellites, there is a possibility of sensor performance degradation (Hu and Le, 2014; Cao et al., 2018). Therefore, regular calibration using the Marine Optical Buoy [MOBY (Clark et al., 1997; Wang et al., 2015)] dataset and validation using the Aerosol Robotic Network - Ocean Color component [AERONET-OC (Zibordi et al., 2006; Mélin et al., 2007)] dataset are essential for these ocean color sensors and have been conducted in several studies (Cannizzaro et al., 2019; Garcia et al., 2020; Pahlevan et al., 2021a). Validations of BPs (Werdell et al., 2009; Cui et al., 2010; Odermatt et al., 2010; Tilstone et al., 2012) and IOPs (Delgado et al., 2019; Fan et al., 2021; Huan et al., 2021) through ocean color imagery have been conducted over regional waters, including ocean, coastal, and inland waters. The consistency of data products from algorithms simultaneously predicting AOPs from multiple satellite missions has been evaluated over different ocean and coastal regions (Cui et al., 2014; Cao et al., 2018; Barnes et al., 2019). Over inland waters, various factors such as the need for atmospheric correction in a more challenging atmosphere due to complex aerosol compositions (compared to oceanic aerosols) and compensation for adjacency effects from nearby lands, produce significantly higher uncertainty for remote sensing (Liu et al., 2024), and data products from multimission datasets have not yet been formally assessed for consistency (Cao et al., 2018; Cao et al., 2023; Pahlevan et al., 2024). A comprehensive assessment that encompasses a wider range of factors, ensuring the reliability and applicability of satellite-derived ocean color data on a global scale, especially for coastal and inland water ecosystems is urgently needed.

The multimission ocean color satellite dataset, along with the extensive hyperspectral field observations from the augmented GLORIA and SeaBASS databases, provides crucial support for analyzing performance of multimission sensor products in both nearshore and inland waters. This study aims at developing a technique based on MDNs for retrieving 10 water quality variables including two BPs and three IOPs at reference bands over coastal and inland waters. The outline of this manuscript is as follows. In Section 2, we discuss the general techniques for obtaining IOPs from ocean color R _rs . Section 3 explains the details of the augmented GLORIA and SeaBASS in situ datasets and the four ocean color satellite sensor datasets used for estimating the 10 BP and IOP variables. Section 4 explains the deep learning method used for the estimation of these water quality parameters from R r s of the four heritage ocean color sensors. Section 5 explains the validation of the developed model with the in situ dataset and multimission ocean color sensor dataset encompassing MODIS/A, MODIS/T, MERIS and VIIRS. In addition to the validation results, the spatial data analysis covers water quality parameters using the near simultaneous measurements of the aforementioned ocean color sensors. Section 6 describes the spatial data analysis. The results are discussed in Section 7, and the study is summarized in Section 8.

In general, R _rs derived from satellite data, can be estimated from the normalized water-leaving radiance (nL _w ) by dividing it by the constant mean extraterrestrial solar irradiance (F ₀ ). This R _rs is associated with other AOPs, which depend on the ambient light field and water constituents, as well as IOPs. Enhancing the relationship between R r s and IOPs in natural waters is fundamental for the accurate estimation of waters constituents from radiometric observations. For this, R _rs can be converted to subsurface remote sensing reflectance r _rs using a relationship (Equation 1) provided as (Lee et al., 2002): r r s λ = R r s λ 0.52 + 1.7 R r s λ (1)

Here, r _rs does not have the effect of the sun being off-zenith associated with it (Werdell et al., 2018). The simplified correspondence between r r s and IOPs revolves around absorption ( a ) and backscattering ( b b ) coefficients, with g representing a complex function dependent on several factors including wavelength, solar angle, wind speed, and water optical properties, given as Equation 2 (Gordon et al., 1988): r r s λ = g λ × b b λ a λ + b b λ (2)

Quantifying the fraction of incident light absorbed within the water column per unit distance, absorption coefficients encompass contributions from pure water a w , phytoplankton a p h , non-algal particles a n a p , and colored dissolved organic matter a c d o m . The combined absorption by non-algal particles (a _nap , also denoted as a _d ) and colored dissolved organic matter ( a c d o m , also known as a _g ) can be referred to as a _dg . The total absorption coefficient (unit: m⁻¹) can be expressed as Equation 3: a λ = a w λ + a p h λ + a c d o m λ + a n a p λ (3)

In clear oceanic conditions (Case 1), a parametrization function can effectively relate Chla to total particulate absorption due to minimal contributions from non-algal particles. However, in coastal and inland waters (Case 2), variations in components like phytoplankton, suspended solids, and colored dissolved organic matter may occur independently from Chla, leading to diverse spectral signatures of R _rs.

The a p h has been related to Chla and phytoplankton biomass (Binding et al., 2008) and can be determined by multiplying their specific absorption coefficients with the concentration of Chla (Roesler and Perry, 1995), as follows: a p h λ = a p h * × C h l a (4)

Here, a p h * in Equation 4 is the phytoplankton-specific absorption coefficient, defined per unit concentration of Chla, that varies widely depending on phytoplankton species composition (Bricaud et al., 1995).

The absorption coefficient for colored dissolved organic matter (see Equation 5) typically follows an exponential decrease as wavelength increases in the near-UV and visible spectral regions, and is observed as (Bricaud et al., 1981): a c d o m λ = a c d o m λ 0 × exp − S × λ − λ 0 (5)

Here, a c d o m ( λ 0 ) denotes the absorption coefficient at the reference wavelength λ 0 (Brezonik et al., 2015). The value of a c d o m ( λ 0 ) is often utilized to characterize the colored dissolved organic matter concentration in a specific water body. Additionally, the spectral slope S (nm⁻¹), which remains independent of λ 0 , signifies the rate at which absorption decreases with increasing wavelength.

The last absorption component, non-algal particles a n a p , generally consists of non-algal organic detritus, living non-algal particles and suspended sediments (Binding et al., 2008). Further, the a n a p λ decreases exponentially with increasing wavelength and the general (Equation 6) is as follows (Roesler et al., 1989; Bricaud et al., 1998; Binding et al., 2008): a n a p λ = a n a p λ 0 × exp − S × λ − λ 0 (6)

Here, a n a p at reference wavelength λ 0 and slope S can be estimated as a function of TSS concentration using an empirical relationship (Babin et al., 2003).

Another IOP is the backscattering coefficient (b _b ), which is typically expressed using two backscattering components (Equation 7): pure water backscattering (b _bw ) and particulate backscattering (b _bp ). The wavelength dependence of b _b is expressed as follows: b b λ = b b w λ + b b p 555 555 λ Y (7)

Here, the values of b _bw are derived from laboratory measurements (Morel, 1974). The b _bp component is influenced by particles in seawater, such as minerals, phytoplankton, and organic matter (Vaillancourt et al., 2004; Woźniak et al., 2019). Y represents the spectral slope of b _bp . Y and b _bp at a reference wavelength of 555 nm can be estimated from R _rs (λ). Most studies use these general (Equations 1–7) to retrieve IOPs from R _rs (Maritorena et al., 2002; Lee et al., 2002; Werdell et al., 2013).

This section elaborates on the detailed methodology utilized in this study for extracting BPs and IOPs from three ocean color missions, as depicted in Figure 5. The approach integrates a MDN to construct a robust deep learning model. To build the model, we employed the augmented GLORIA dataset integrated with IOPs for training and testing purposes. For model evaluation, we utilized a matchup dataset sourced from three different sensors: MODIS, MERIS, and VIIRS, aligned with corresponding SeaBASS datasets. This section presents a detailed overview of the MDN model development and the procedures for processing satellite images.

The MDN architecture achieves low median retrieval errors (<45%) for all 10 variables of BPs and IOPs from in situ R _rs . The test results on 20% of the in situ augmented GLORIA samples are shown in the Appendix (refer to Supplementary Appendix Figures SA1-A3, and performance metrics in Supplementary Appendix Tables SB3-B6) at the spectral resolutions of the MODIS, MERIS, and VIIRS sensors. This prior analysis facilitates the evaluation of the model’s performance with satellite datasets through matchup analysis and the results of the leave-one-out cross-validation are discussed in the following sub-sections.

In the previous section, we found a reasonable alignment between MODIS, MERIS, and VIIRS products and the MDN-based BPs and IOPs, particularly when compared to in situ datasets across coastal and inland waters. This section elucidates the efficacy of utilizing spatial data to comprehend the dynamic fluctuations in geophysical products within coastal and inland water bodies. Specifically, we focus on the Chesapeake Bay to analyze the MDN-based spatial products derived from MODIS, MERIS, and VIIRS. Situated along the East Coast of the United States, Chesapeake Bay encompasses highly productive waters, exhibiting a spectrum of phenomena including in-water blooms, suspended sediments, and dissolved materials across the bay and coastal regions.

Figure 11 presents MODIS/T (a) and MERIS (b) ocean color products showcasing the spatial distribution of BPs and IOPs on 27 December 2005, post-application of the MDN method. Leveraging the MDN retrievals, these images facilitate a clearer understanding of the spatial distribution of water quality parameters. Regions depicted in red indicate a high dominance of the corresponding water quality parameter, while those in blue signify a lower one. Figures 11a1, 11b1 illustrate the spatial variation of MDN-derived Chla as captured by MODIS/T and MERIS, respectively. The MODIS/T image vividly highlights the dominance of Chla near the shoreline in the northern and western parts of the bay. Conversely, the MERIS image exhibits less variation compared to MODIS, with less discernible features. The spatial variation of MDN-derived TSS, is depicted in Figures 11a2, 11b2 for MODIS/T and MERIS, respectively. In contrast to Chla, the dynamic changes in TSS distribution are readily discernible across the entire bay area and coastal waters in both MODIS/T and MERIS images.

Similar to the BPs, the spatial distribution of the IOPs also distinctly reveals their variation in these optically complex waters. Figures 11a3, 11b3 highlight the presence of a _ph at 443 nm near the shoreline of the bay. The dominant features of the a _ph distribution are clearly depicted in the MODIS/T images, whereas they are less evident in the MERIS images. In Figures 11a4, 11b4, the distribution of a _cdom at 443 nm is distinctly visible, demonstrating that both MODIS/T and MERIS effectively capture the spatial distribution of dissolved organic matter over the bay area. This organic matter dissolves into the water due to tidal movements between the bay and the sea. The MODIS/T-acquired a _cdom image over the bay area aligns consistently with the MERIS-acquired a _cdom image. However, notable discrepancies arise offshore, where the two satellite images display significant variation. The spatial distribution of another crucial IOP variable, a _nap at 443 nm, is depicted in Figures 11a5, 11b5. Remarkably, the a _nap products obtained from both MODIS/T and MERIS sensors exhibit close consistency. Elevated levels of a _nap are observed predominantly along the shoreline of the bay and in coastal waters. This analysis distinctly elucidates that the BPs and IOPs derived from the MODIS/T and MERIS sensors exhibit consistency for TSS and a _nap , while slight variations are observed for other BPs and IOPs.

Figure 12 depicts MODIS/A (a) and VIIRS (b) ocean color products highlighting the spatial distribution of BPs and IOPs over Chesapeake Bay on 26 December 2018. Surprisingly, all ocean color products examined exhibit matching spatial distributions for both sensors. Evaluating the MDN-derived BPs Chla and TSS, as well as the IOPs, for consistency and spatial variability, clear patterns can be discerned. Figures 12a1, 12b1 illustrate the MDN-derived Chla levels detected by MODIS/A and VIIRS, respectively, denoting a notable concentration in the bay area and particularly along its shoreline. Conversely, Chla concentrations appear considerably lower in coastal and offshore waters. Noteworthy peaks in Chla are observed exclusively along the bay’s shoreline. The spatial distribution of MDN-derived TSS, is depicted in Figures 12a2, 12b2. Both MODIS/A and VIIRS ocean color products exhibit heightened TSS levels along the bay’s shoreline, with a consistent pattern, albeit slightly lower magnitudes in VIIRS compared to MODIS/A.

MDN-derived IOPs, such as a _ph , are coherent across both MODIS/A and VIIRS ocean color products, as demonstrated in Figures 12a3, 12b3. These images highlight an abundance of phytoplankton pigments within bay waters, especially in proximity to the shoreline, with minimal presence in offshore regions. This trend is replicated across both MODIS/A and VIIRS datasets. Similarly, Figures 12a4, 12b4 showcase the dominance of a _cdom in bay waters, with negligible quantities observed offshore. The spatial distribution of a _cdom , as estimated by the MDN, remains largely consistent across both platforms. Contrary to other IOPs, a _nap exhibits an uneven distribution within the bay zone, as depicted in Figures 12a5, 12b5. Non-algal particles are predominantly concentrated near the bay’s shoreline, coinciding with areas of heightened TSS levels (as shown in Figures 12a2, 12b2). Notably, the spatial patterns of MDN-derived a _nap remain highly similar for both MODIS/A and VIIRS imagery.

The consistency in spatial patterns between MODIS/A and VIIRS images suggests reliable performance in capturing ocean color variability over Chesapeake Bay. While minor differences exist in magnitude, overall trends remain consistent between the two sensors. Notably, both MODIS/A and VIIRS images show similar spatial distributions for MDN-derived BPs and IOPs. The comparative analysis indicates that MODIS/A and VIIRS sensors provide consistent and reliable ocean color products for monitoring the Chesapeake Bay, even though slight discrepancies exist, particularly with MERIS images.

Recent studies have revealed the immense potential of satellite-based hyperspectral observations in enhancing our understanding of aquatic ecosystems on a global scale (O’Shea et al., 2023). However, the reliability of product retrievals tends to diminish with a multispectral sensor or equivalent reduction in the number of bands utilized. In our investigation, we scrutinized three different multispectral sensors - MODIS, MERIS, and VIIRS - each characterized by a distinct number of multispectral bands. Surprisingly, despite this variance, the accuracy of the MDN-estimated BPs and IOPs exhibited minimal discrepancy across the sensors, as depicted in the training results (first column, Figure 13), leveraging a substantial subset of 20% of the augmented GLORIA in situ dataset. Notably, we demonstrated that achieving a 20%–25% error margin in all BPs and IOPs utilizing MDNs is possible on in situ measured data. Slight disparities observed in the accuracy of BPs and IOP variables among the sensors can be attributed to differences in the available multispectral bands within each sensor.

The heritage AC model represents a significant advancement in the field, delivering improved R r s products compared to prior studies (Pahlevan et al., 2024). However, challenges persist, particularly in highly productive aquatic environments where negative retrievals are frequently encountered. These discrepancies introduce uncertainties, especially in extrapolating aerosol contributions to the blue region of the R r s spectrum (Frouin et al., 2019). Such uncertainties reverberate throughout downstream products like BPs and IOPs, as evidenced in this study. For example, our analysis reveals that residual errors in R r s lead to approximately 50% overestimation across key parameters such as Chla, TSS, a _cdom , a _nap , and a _ph (second column, Figure 13) for the same GLORIA dataset. Here, satellite R r s is used as the input for the MDN model. The geographical distribution of GLORIA matchups (+/−3 h) with their year wise distribution for these ocean color missions are shown in the appendix (Supplementary Appendix Figures SA6, A7). The scatterplots and performance metrics for GLORIA matchups with MODIS, MERIS, and VIIRS are presented in the appendix as well (Supplementary Appendix Figures SA8, A11 and Supplementary Appendix Tables SB6-B9) and illustrate how issues in atmospheric correction directly impact R r s and consequently influence the water quality variables.

The results highlight that MODIS use of both NIR and SWIR bands enables effective atmospheric correction, minimizing atmospheric effects. While VIIRS employs a similar approach, its limited visible bands restrict its ability to resolve complex optical properties in aquatic environments, affecting the retrieval of inherent optical properties (IOPs) like a _cdom and a _nap . MERIS, with a broader visible spectrum, offers better spectral resolution, but relies solely on NIR bands for atmospheric correction, potentially introducing biases under certain conditions. These limitations in VIIRS and MERIS atmospheric correction and band configurations may lead to lower a _cdom and a _nap estimates. The validation dataset for VIIRS, particularly for a _cdom and a _nap , is limited, potentially affecting the robustness of these specific validations. To address this, the Appendix includes MDN performance evaluations for the GLORIA dataset and VIIRS sensor, offering initial insights into model capabilities. While matchups for a _cdom and a _nap are sparse, the datasets cover diverse global waters, supporting the MDN’s applicability across multiple variables. Additionally, we highlight the model’s performance on other water quality indicators, such as Chla, TSS, and a _ph , where data coverage is more comprehensive, providing a broader and more reliable evaluation of its capabilities.

The third column (Figure 13) shows the performance metrics for the SeaBASS matchups with the three satellite sensors. As we discussed in Section 5.1, the error still increased approximately 75%. It is plausible that these inaccuracies stem not only from the model and atmospheric correction, but also from factors such as inadequate vicarious calibrations, residual biases in TOA measurements (Ibrahim et al., 2018), and adjacency effects (Sterckx et al., 2011). Our findings strongly suggest that the MDN is particularly susceptible to errors in atmospheric correction (Zolfaghari et al., 2023) and adjacency effects. Addressing these challenges requires a concerted effort of both refining the models and emphasizing the critical need for accurate atmospheric correction in aquatic remote sensing studies.

This study introduces an MDN-based inversion method tailored for the concurrent estimation of Chla, TSS, a _cdom , a _nap , and a _ph from R r s across a range of multispectral ocean color sensors, applicable to coastal and inland water bodies. By training the MDN model using additional IOP data aligned with 80% of the augmented GLORIA dataset, better accuracy is achieved. Notably, the model demonstrates enhanced or comparable accuracy (ranging from 15% to 40% across all products) when validated against a larger and more diverse in situ dataset compared to prior dedicated MDN approaches. Validation of the SeaBASS dataset against MODIS, MERIS, and VIIRS satellite datasets through matchup analysis further underscores the robustness of the model, with accuracy ranging from 40% to 70% across all BPs and IOPs.

Leveraging a regional leave-one-out cross-validation approach, the MDN exhibits superior generalization performance across various sensor platforms and water types. Furthermore, this study elucidates the suitability of MDN models for specific regional water bodies and identifies sensors that offer better accuracy for individual BPs and IOPs. The results of the leave-one-out analysis reveal minimal discrepancies among the retrievals from the three ocean color sensors. When applied to satellite images, slight differences between MODIS, MERIS, and VIIRS are observed, primarily attributed to atmospheric correction algorithm nuances and sensor band disparities.

Addressing the substantial uncertainty stemming from atmospheric correction underscores the potential for significant accuracy enhancements in IOP retrieval through improved remote sensing reflectance retrieval techniques. Future refinements to the MDN approach could incorporate additional environmental and physical variables as input features, facilitating phytoplankton type/species discrimination within specific geographic regions. In this study, the MDN approach was applied exclusively to heritage ocean color sensors, driven by the absence of a matchup validation dataset for recent ocean sensors. The availability of in situ data presents an opportunity to expand this methodology to recent multispectral sensors such as VIIRS onboard NOAA-20 and NOAA-21 and OLCI onboard Sentinel-3A and Sentinel-3B. These advancements hold promise for enhancing monitoring capabilities of water bodies.