Simultaneous retrieval of 13 phytoplankton pigments from Sentinel-3 OLCI using convolutional neural networks and transfer learning

1 Introduction

Phytoplankton pigment concentrations are generally measured through in situ water sampling and subsequent analysis via High Performance Liquid Chromatography (HPLC), a pressure-driven liquid chromatographic technique in which compounds are separated according to their differential retention between a liquid mobile phase and a particle-packed stationary phase, allowing reliable identification and quantification of non-volatile compounds. HPLC permits the contemporaneous retrieval of the concentration of numerous pigments with an average percentage error, estimated as the relative difference between analyses of the same sample at different laboratories, on the order of 20% (Canuti, 2023). However, both the sampling activities and the analyses are extremely time-consuming and expensive, which limits the spatio-temporal characterization of large basins. Alternative to in situ sampling, remote-sensing can provide frequent observations with high spatial coverage and resolution. A number of algorithms have been developed in the recent years to retrieve pigment concentrations from remote-sensing measurements based on various techniques. Among those, Pan et al. (2010) applied empirical algorithms to multi-spectral observations provided by the Moderate Resolution Imaging Spectroradiometer (MODIS) or the Sea-viewing Wide Field-of-view Sensor (SeaWiFS) in the northeastern coast of the United States. Differently, Bracher et al. (2015) developed a set of empirical orthogonal functions and applied a multiple linear regression based on co-located Medium Resolution Imaging Spectrometer (MERIS) observations in the Atlantic Ocean, while Moisan et al. (2017) performed a linear matrix inversion based on the phytoplankton absorption spectra, in turns retrieved from multi-spectral observations provided by MODIS. Chase et al. (2017) decomposed hyper-spectral absorption spectra, in turns derived by hyper-spectral radiometric measurements, as a sum of Gaussian functions and used their inverted amplitude to estimate pigment concentrations. Noteworthy, El Hourany et al. (2019) applied self-organizing maps to a satellite-based dataset including multi-spectral radiometry, chlorophyll-a and sea surface temperature. Sun et al. (2022) employed spectral decomposition relying on multi-spectral SeaWiFS and MODIS observations. Kramer et al. (2022) focused on the residual of hyper-spectral remote sensing reflectance measurements after removing the principal features associated to chlorophyll-a, dissolved and non-pigmented matter, and particulate matter. Recently, Li et al. (2023) applied a deep-learning model to satellite measurements provided by MODIS.

All these studies provide alternative methods to retrieve pigment concentrations from space. However, from the radiometric point of view, most of them, with the exception of Chase et al. (2017) and Kramer et al. (2022), rely on satellite matchups. These are characterized by a relatively limited radiometric accuracy and high sensitivity to spatial inhomogeneities, which may eventually result in an increased uncertainty of the retrieved pigment concentrations. Moreover, concerning the spectral sampling, none exploits the improved spectral features, i.e., the specific additional bands, of the Ocean and Land Color Instrument [OLCI, Donlon et al. (2012)] which, with respect to those of previous sensors such as MODIS, MERIS, or SeaWiFS, may ease the pigment retrieval process. Within this context, the main objective of this study is fourfold. First, it provides the first systematic and controlled benchmark of five machine learning architectures (Random Forest, eXtreme Gradient Boosting, Dense Neural Network, one-dimensional Convolutional Neural Network, and Bidirectional Long Short-Term Memory Network) for the simultaneous retrieval of 13 diagnostic pigments from OLCI-equivalent multispectral radiometry, leveraging the sensor’s specific band configuration, which has remained largely unexplored in previous studies. Second, it introduces and validates a novel two-phase deep learning training strategy designed to prevent spectrally dominant pigments from compromising the retrieval of less abundant, or less represented, accessory pigments within a joint multi-output framework. Third, it explicitly quantifies the value of high-quality regional in situ radiometry as a pre-training foundation for global satellite pigment retrieval by comparing a fine-tuned model with an identical architecture trained from scratch using satellite matchup data. Fourth, it develops and validates a spectrally reduced five-wavelength version of the best-performing model, enabling its application to long-term merged satellite records such as Copernicus-GlobColour (Maritorena et al., 2010) or the ESA Ocean Colour Climate Change Initiative (Sathyendranath et al., 2023).

2 Materials and methods

2.1 In situ dataset

This study relies on a set of coincident measurements of remote sensing reflectance ( R r s ) and concentrations of a set of pigments of interest to the ocean color community carried out between 2014 and 2021 in European waters (185 measurement stations) in the framework of the Bio-Optical mapping of Marine optical Properties [BiOMaP, Zibordi et al. (2011); Zibordi and Berthon (2024)]. The geographic distribution of the in situ measurements included in the study is shown in Figure 1.

Radiometric measurements were performed using TriOS (Rastede, DE) RAMSES hyper-spectral radiometers operated in either above-water or in-water single-depth configurations, or alternatively, BioSpherical Instruments (BSI, San Diego, CA) Compact-Optical Profiling System (C-OPS) multi-spectral radiometers operated in in-water free-fall profile configuration. In all cases, data acquisition strictly followed the protocols detailed in Zibordi et al. (2019) to minimize the radiometric uncertainty. TriOS RAMSES operate in the nominal spectral range of 320–950 nm, with spectral sampling and resolution of approximately 3.3 and 10 nm. Data were calibrated and processed as detailed in Talone et al. (2016), Talone et al. (2018), Talone et al. (2024) to produce R r s spectra with a regular spectral sampling of 2 nm in the 400–800 nm spectral interval. For each of the spectra, OLCI-like multi-spectral observations were obtained through linear interpolation of the hyper-spectral measurements. On the contrary, BSI C-OPS incorporates multi-spectral radiometers collecting simultaneous measurements at 19 center-wavelengths. Among those, 11 C-OPS center-wavelengths were selected to exactly match OLCI’s one, so no spectral interpolation was needed in this case. Measurements were performed and processed as described in Zibordi et al. (2011). Basic statistics of the in situ radiometric dataset such as average value ( m ), standard deviation ( σ ), maximum ( m a x ) and minimum ( m i n ) values are summarized in Table 1 for each of the selected central-wavelengths ( λ ).

TABLE 1. Summary statistics of the radiometric measurements included in the in situ dataset at the various wavelengths.

Complementing Table 1, spectral plots of the R r s measurements are shown in Figure 2. Similar to Talone et al. (2018), blue, green, and orange solid lines indicate spectra corresponding to Case-1 waters, Case-2 waters likely dominated by moderate-to-high suspended sediment concentrations (Case-2a), and Case-2 waters likely dominated by high CDOM (Chromophoric Dissolved Organic Matter) concentrations (Case-2b), respectively. R r s values are always lower than 0.015 s r − 1, and generally lower than 0.01 s r − 1, indicating mainly mesotrophic waters with medium-to-low concentrations of sediments; spectral maxima are relatively evenly distributed, suggesting the presence of spectra from a variety of water types, from the quasi-oligotrophic waters of the Balearic Sea to the complex waters characterized by moderate presence of sediments and/or dissolved organic matter of the Adriatic Sea and the Black Sea.

Additionally, for each of the center-wavelengths, the normalized histogram of R r s is shown in Supplementary Figure S1 together with the corresponding Kernel Density Estimation (KDE) over-plotted in red. Distributions are generally multi-modal, reflecting the diversity of water types observed in Figure 2, and mostly asymmetrical, suggesting the need for a log-transformation of the data before the algorithms training phase.

Regarding pigment concentrations, water samples were collected coincident with radiometric observations, filtered using 0.7 μ m GF/F filters, and successively analyzed by HPLC. The result of the HPLC analysis, performed by the Danish company DHI (Hørsholm, DEN) in accordance with Hooker et al. (2012), was the quantification of the concentration ( C ) of 13 pigments, namely Chlorophyllide a (chlide), Chlorophyll a (chla), Chlorophyll b (chlb), Chlorophyll c1 + c2 (chlc12), Fucoxanthin (fuco), 19′-Hexanoyloxyfucoxanthin (hex), 19′-Butanoyloxyfucoxanthin (but), diadxanthin (diad), Alloxanthin (allo), Diatoxanthin (diato), Zeaxanthin (zea), Beta-carotene (caro), and Peridinin (peri). Similarly to the radiometric measurements, average value ( m ), standard deviation ( σ ), and minimum ( m i n ) and maximum ( m a x ) values for each of the pigments included in the input dataset are shown in Table 2, whereas normalized histograms and corresponding KDEs are displayed in Supplementary Figure S2, with blue bars and red lines, respectively.

TABLE 2. Summary statistics of the pigment concentrations ( C ) included in the in situ dataset.

Pigments exhibit widely different distributions, with average values and standard deviations spanning two orders of magnitude. Namely, average values are comprised between 0.017 and 0.811 m g m − 3 whereas standard deviations range between 0.017 and 1.674 m g m − 3, indicating different levels of retrieval complexity depending on the pigment. Similarly to R r s measurements, KDE plots exhibit asymmetrical distributions for most of the pigments, reinforcing the need of a log-transformation of the input data. For both pigments and radiometric datasets, a lower bound equal to the 5 t h percentile of each variable is imposed prior to log-transformation. Values below this threshold are replaced with the corresponding percentile value, preventing numerical instabilities associated with computing the logarithm of near-zero values.

2.2 Satellite match-ups dataset

Complementing the in situ observations, an independent dataset of co-located pigment observations and satellite-derived R r s (matchups) with global coverage was used to extend the applicability of the algorithm. The pigment concentrations dataset was assembled from multiple open-access in situ databases, including the MARine Ecosystem DATa [MAREDAT, Peloquin et al. (2013)], the NASA bio-Optical Marine Algorithm Dataset [NOMAD, Werdell and Bailey (2005)], and the TARA dataset (Karsenti et al., 2011; SeaBASS, 2009), while the Global Ocean (Copernicus-GlobColour) product provided by the Copernicus Marine Service (European Union-Copernicus Marine Service, 2022) was selected as source of the satellite dataset. This latter choice was mainly motivated by (i.) the availability in the same product of both the OLCI-like (11 center-wavelength) and SeaWiFS-like (5 center-wavelengths) spectral configurations (ii.) its global coverage that maximizes the number of match-ups permitting a more robust performance analysis, and (iii.) the availability of a documented Application Programming Interface (API) for its download. It is noted that although the OC-CCI dataset Sathyendranath et al. (2023) could also have been used for this purpose, it was not explicitly evaluated in this study. Nevertheless, due to its substantial spectral overlap with GlobColour, no compatibility issues are expected, and the implementation code, publicly available, allows straightforward application to OC-CCI product.

Two datasets were defined based on both the 11- and 5-center-wavelength products, respectively. These were built considering a spatial radial window of 0.06 × 0.06 degrees centered on the location of the HPLC sample and a temporal window of ± 1 day. A total of 668 and 7,287 matchups were identified for the 11- and 5-center-wavelength satellite-based datasets, with a geographical distributions illustrated in Figure 3.

Spectral plot, histograms, and main statistics of both satellite-based datasets, shown in Supplementary Figures S5–S8 and Supplementary Tables S1–S4, differ notably from those of the in situ dataset. For instance, maximum values observed in the in situ radiometry never exceed the value 0.015 s r − 1, while the maximum values found in both the 11- and the 5-center-wavelength datasets are consistently above 0.019 s r − 1. Significant differences are also noted in pigments statistics, with both average values and standard deviations showing very diverse values suggesting the inclusion of waters exhibiting a wider range of optical properties. Additionally, some pigments are not equally represented in the satellite-based datasets compared to the in situ dataset, with the most extreme case being Beta-carotene caro, which is entirely absent from the satellite-based data of 11-wavelengths.

2.3 Model description

Two categories of models were evaluated: classical Machine Learning and Deep Learning models. As summarized in Table 3, the first category includes Random Forest (RF) and Extreme Gradient Boosting (XGB), which are among the most widely exploited models in the marine remote sensing community (e.g., Medeiros et al., 2015; Cheng et al., 2021; Xiao et al., 2022; Chen et al., 2024; Kim et al., 2025). Random Forest is based on the generation and averaging of a series of random decision trees; its strengths include robustness, high performance in scenarios with limited input datasets, and a straightforward interpretability. Also based on the generation of decision trees, Extreme Gradient Boosting differs from RF by building decision trees sequentially, where each new tree corrects the errors of the previous ones. Beyond the advantages shared with RF, key features of XGB include its speed, L1 and L2 regularization techniques, and its capability to handle missing values. Regarding Deep Learning, the tested models include a Dense Neural Network (DNN), a one-dimensional Convolutional Neural Network (CNN), and a more complex Bidirectional Long Short-Term Memory Neural Network (BiLSTM). The first consists of a concatenation of dense layers (fully connected) with dropout layers in between. The activation function of the dense layers is a Rectified Linear Unit (ReLU), except for the last one, which applies a standard linear activation. The CNN is built on one-dimensional convolutional layers, a structure that provides robust generalization capabilities and mitigates overfitting when processing sequentially correlated input data. Finally, the BiLSTM model is a recurrent neural network consisting of several bidirectional LSTM layers, typically used in natural language processing and problems with a notable sequential correlation in the input data features.

TABLE 3. Models evaluated and their classification into classic machine-learning and deep-learning models.

2.4 Pre-processing, training, and testing

A series of transformations were applied to the dataset to facilitate the application of ML techniques. Firstly, to improve models’ training capabilities, both pigment concentrations and radiometric measurements were log-transformed, resulting in an increased normality and more stable variance (Campbell, 1995; O’Reilly et al., 1998). Additionally, a stratified data-splitting approach with a 15% factor was employed for the training and test datasets. The retrieval of diagnostic pigment concentrations, in fact, is a multi-output problem, which entails the need to stratify the output into different classes, or clusters, prior to data splitting to ensure that the data distribution is preserved across the various splits. To achieve this, K-means clustering was applied to the model outputs prior to data splitting. K-means partitions the data into K non-overlapping clusters by minimizing within-cluster variance, ensuring that spectrally similar samples are grouped together. Based on the number of observations (185), the number of clusters (K = 15) was chosen as the maximum value ensuring that each cluster contained at least two observations, a requirement for subsequent stratified splitting. Lastly, to overcome the risk of overfitting due to the limited size of the training dataset and with the objective of introducing realistic variability without altering the underlying data distributions, two data augmentation techniques were applied, namely.

• Sample mixing: Pairs of samples were linearly mixed with a 90%–10% weighted average to simulate intermediate radiometry and pigment combinations; and
• Noise injection: Random noise sampled from normal distributions of mean 0 and standard deviations 0.05 and 0.17 was applied to simulate a 5% and 17% measurement uncertainties of radiometric data and pigment concentration. These figures are consistent with the results of previously published field and laboratory uncertainty assessment studies, namely the 5% radiometric uncertainty target in the blue-green spectral region is widely accepted in the Satellite Ocean Color Validation Community (e.g., Zibordi et al., 2014), while the 17% uncertainty assigned to HPLC data results from the average value of the absolute percent differences between pigments concentrations obtained in a series inter-laboratory round-robin exercises described in Canuti (2023).

These augmentations expanded the training set by a factor of 100. Examples of the augmented spectra are shown in Supplementary Figure S9. The synthetic spectra are represented by thin gray lines, while the thick blue, green, and orange lines correspond to representative measurements from the Case-1, Case-2a, and Case-2b in situ water datasets, respectively. Based on the augmented datasets, the best hyper-parameters for each model were selected through cross-validation, and the predictive capacity of the models was evaluated against the test set to provide estimates of generalization errors. To increase the robustness of this estimation, the entire procedure was repeated six times, including stratified splitting, data augmentation, training, and testing. This step is common to both ML and DL models. A schematic of the pre-processing, training, and testing methodology is shown in Figure 4.

2.5 Metrics

The performance of the predictive models was assessed based on a set of consolidated metrics, specifically:

• • The coefficient of determination (R ², Equation 1) and the Mean Squared Error (MSE, Equation 2), computed in the logarithmic space of the pigments as:

R 2 = 1 − ∑ i n log 10 C i ′ − log 10 C i 2 ∑ i n log 10 C i − m log 10 C 2 (1)

M S E = 1 n ∑ i n log 10 C i ′ − log 10 C i 2 , (2)

where subscript i identifies each of the n sample points, C and C ′ indicate measured and predicted pigment concentrations, and m is the average value. These provide information regarding the optimization of the model’s hyper-parameters, which, as mentioned, is performed in the logarithmic space. Notably, the score applied in both the training and the cross-validation of the model is based on the M S E.

• Complementary, the Mean Absolute Error (MAE, Equation 3) and the Mean Absolute Percentage Error (MAPE, Equation 4) were computed in the actual pigment space to quantify the dispersion of the predicted values around the measured ones in absolute and percentage terms:

M A E = 1 n ∑ i n | C i ′ − C i | , (3)

M A P E = 100 n ∑ i n | C i ′ − C i | C i (4)

• The Mean Error (ME, Equation 5) and the Mean Percentage Error (MPE, Equation 6) were also computed in the actual pigment space to quantify the systematic bias of the model, expressed in absolute and relative terms, respectively:

M E = 1 n ∑ i n C i ′ − C i , (5)

M P E = 100 n ∑ i n C i ′ − C i C i (6)

As noted, satellite-based matchup dataset combines HPLC measurements from multiple publicly available datasets with heterogeneous pigment suites, resulting in missing values for pigments not reported in certain databases. To accommodate this, sample observations with unknown values ( N a N ) are excluded from the computation of all metrics, enabling both training and evaluation on incomplete datasets.

2.6 Deep learning architecture

Deep Learning models may be more sensitive to overfitting, as their capacity to train and synchronize large numbers of parameters may increase the risk of capturing undesired patterns specific to the training dataset. Moreover, pigments exhibited highly variable magnitudes and distributions (Table 2; Supplementary Figure S2), indicating diverse levels of retrieval complexity. While some pigments may overfit early during training, others could benefit from longer training with a higher number of epochs and more complex architectures. To mitigate the risk that optimizing performance for certain pigments might negatively affect others, the DL models employed a two-phase training strategy consisting of an initial pigment-independent training phase, followed by a simultaneous multi-pigment phase. This allows for the optimization of training for each pigment while simultaneously enabling the model to leverage inter-output dependencies, information that remains inaccessible to models such as RF or XGB. A schematic of the models in the first training phase is shown in Figure 5. Each module was trained to infer a single pigment. Depending on the main layers exploited, three different models appear in this phase: a DNN model with dense layers, a CNN model with convolutional layers, and a BiLSTM model with Bidirectional LSTM layers. All of them include an initial normalization layer to standardize the input data according to the training set. A complete 3-fold cross-validation process was performed to find the best hyper-parameters across previously designed model architectures (see Supplementary Tables S5–S12, for a detailed architectures definition). An early-stopping mechanism with a patience of two was employed to prevent overfitting.

Following the first single-pigment training phase, a second model was developed based on the concatenation of the modules from Phase 1 and incorporating a dense head layer at the output. A schematic of the second training phase model is shown in Figure 6. Training was carried out for 300 epochs, with the weights obtained during Phase 1 kept frozen and only the head layer set as trainable. This was followed by an additional 100 epochs during which all layers were unfrozen and optimized jointly with a lower learning rate. By forcing the model to predict all pigments from the fixed modules, this step permit to obtain a synchronized and interdependent response.

3 Results and discussion

The performance of all evaluated models is presented in this section on a pigment-by-pigment basis. It is noted that all the results reported in this section correspond to the final experiments conducted on the held-out test sets.

3.1 Comparative analysis of the ML models

Performance metrics of all ML models are summarized in Table 4 with the best results highlighted in bold.

TABLE 4. Metric values (R 2, M S E, M A E, M A P E, M E, and M P E) on a pigment-by-pigment basis for each of the tested ML models. Best metrics are highlighted in bold.

The results show a clear performance hierarchy in the R 2 metric, with CNN emerging as the dominant architecture, while BiLSTM, RF, and XGB remain competitive for selected pigments and DNN generally shows lower performance. A strong coherence is obtained between R 2 and M S E across almost all pigments. In 10 out of 13 cases, CNN has the highest R 2 and also achieves the lowest M S E. This consistency confirms that the models are successfully minimizing the loss function during training. Retrieval performance varies substantially among pigments, likely due to their specific spectral features as well as their abundances in the training dataset. Three distinct groups can be identified based on R 2. The first group characterized by R 2 larger than 0.85 includes c h l a, c h l c 12, f u c o, d i a d, c a r o, and p e r i. Within this group, c h l a and f u c o benefit from strong and spectrally distinctive absorption features well covered by OLCI band, particularly around 442, 490, and 673 nm, combined with the widest concentrations ranges in the input dataset, providing the model with robust statistical leverage. The good performance of c h l c 12, d i a d, and c a r o, is instead likely explained by a combination of absorption features that are partially captured by the OLCI bands and their stoichiometric co-variation with c h l a and f u c o, whose well-constrained retrieval likely propagates predictive skill to these co-varying pigments. Lastly, despite substantial spectral overlap between its absorption features and those of f u c o, p e r i also achieves a high R 2, possibly aided by its wide concentration range in the input dataset. In general, for those pigments, both XGB and CNN consistently achieve very high R 2 and M A P E on the order of the absolute percent differences obtained in published inter-laboratory HPLC round robin exercises (Canuti, 2023), suggesting that their retrieval has reached a maturity level potentially suitable for operational use. The second group, with values of R 2 comprised between 0.6 and 0.85, includes c h l i d e, c h l b, h e x, a l l o, and d i a t o. These pigments are characterized by absorption features that in several cases strongly overlap with those of first-group pigments (e.g., c h l i d e versus c h l a, or h e x versus f u c o) and are generally both less represented and characterized by lower concentrations in the input dataset. For this group, performance is more model-dependent, with deep learning approaches often outperforming classical approaches, indicating that learned spectral feature extraction, likely exploiting spectral shape/gradients, offers an advantage over tree-based methods in case of weaker spectral signatures. The third group with R 2 lower than 0.60 comprises b u t and z e a. In this latter case less accurate results are obtained independently from the model applied. This specific behavior is possibly due to the combination of relatively poor representation in the input dataset, low average pigment concentrations, and high spectral overlap with dominant pigments. Although specific values of the various coefficients of determination vary among different studies, similar groups were obtained in Pan et al. (2010), Bracher et al. (2015), Moisan et al. (2017), El Hourany et al. (2019), and Stock and Subramaniam (2020), indicating inherent difficulties in the retrieval of certain pigments, especially b u t and z e a [except in El Hourany et al. (2019)]. An additional confirmation is given by the analysis of M A P E and M P E with respect to M A E and M E. Pigments such as a l l o, d i a t o, z e a and p e r i consistently exhibit worse M A P E and M P E but concurrently show very low M A E and M E values that approach analytical detection limits. The M E is predominantly negative across all pigments and models. This indicates a systematic underestimation bias in the ML models. Opposite sign is observed for M P E. This seemingly contradictory result implies a skewed error distribution. The models likely slightly underestimate the abundant samples (driving negative M E) but significantly overestimate the very low concentration samples (driving positive M P E). The CNN model demonstrates superior performance across most of the error and variance metrics. Specifically, CNN accounts for the highest coefficients of determination and the lowest M S E for 10 of the 13 target pigments. CNN also exhibits the lowest average M A P E, reaching 54.6%, while its average M P E remains competitive at 26.0%. Concerning M E, i.e., an estimate of the systematic bias, XGB shows the best performance with an average absolute bias of approximately 0.013 m g m − 3. Despite its superior complexity, BiLSTM is not able to outperform CNN overall, although it reaches the best R 2 and M S E for h e x and a l l o and remains competitive for several other metrics. The strong performance of the CNN, which is specifically designed to capture spatially localized features, e.g., affecting measurements performed at a short number of consecutive center-wavelengths, indicates that the retrieval of pigment concentration is primarily driven by local spectral information. In particular, distinct sub-patterns within specific wavelength ranges appear to be more informative for pigment dynamics than global or long-range spectral dependencies. The model’s ability to extract these localized features across the spectrum is more effective than the long-range temporal/spectral dependencies modeled by BiLSTM, or the simple feature combination of RF/XGB. The comparison between RF and XGB, both prominent tree-based ensemble methods, reveals a classic trade-off between overall accuracy and systematic bias in the aggregated pigment retrieval results. Among the two tree-based models, RF emerged as marginally superior in terms of average R 2 and M S E, suggesting that its bagging approach is effective at reducing prediction variance and improving robustness. Conversely, XGB exhibited the lowest systematic bias ( M E ) across all tested models. This makes XGB the preferred choice for applications requiring accurate, large-scale estimation with minimal bias, while RF remains a robust baseline for maximizing tree-based prediction accuracy. Overall, CNN is the recommended architecture for maximizing retrieval accuracy ( R 2 ) across the widest range of pigments. However, XGB remains a robust alternative, particularly valuable for its lower systematic bias and computational efficiency on major pigments. A visual synthesis of the CNN pigment-by-pigment performance is shown in Figure 7, where each bubble represent a pigment, horizontal and vertical axes indicate R 2 and the M P E, and the bubble color encodes the M A P E. This visualization complements Table 4 by jointly displaying explained variance, systematic relative bias, and relative dispersion for the best-performing architecture. Overall performance generally decrease from the bottom-right to the top-left corner of the figure.

The bubble representation clearly confirms the pigment-dependent retrieval hierarchy already observed in Table 4. The best-performing pigments, including c h l a, c h l c 12, f u c o, d i a d, and c a r o, appear toward the high-R 2 region and generally exhibit moderate M P E and M A P E. Conversely, b u t and z e a remain in the lower-R 2 region, confirming their higher retrieval difficulty. The case of a l l o is noteworthy: despite an intermediate R 2, it shows very large relative errors, reflecting the sensitivity of percentage-based metrics to very low pigment concentrations.

3.2 Extension to satellite observations

Satellite observations introduce a number of additional sources of uncertainty with respect to in situ ones, primarily related to atmospheric correction, adjacency effects, and sub-pixel heterogeneity. The atmosphere often accounts for more than 90% of the signal measured by a satellite sensor, so that any uncertainty in its characterization largely affects the estimation of the bottom-of-the-atmosphere remote sensing reflectance, and, consequently, of the pigment concentrations derived from it [Mélin (2019)]. The proximity of land or clouds can further perturb the radiance collected by the sensor through the so-called adjacency effect: depending on the magnitude and spectral characteristics of the surrounding land or cloud reflectance, this contamination can contribute more than 100% of the unperturbed water-leaving reflectance (Bulgarelli and Zibordi, 2018). Additionally, sub-pixel spatial heterogeneity may further increase the measurement uncertainty by introducing a representativeness mismatch between spatially integrated satellite observations and point-scale in situ measurements (Pahlevan et al., 2016). This is particularly pronounced in dynamically active regions such as upwelling fronts, river plumes, and eddies, where sharp optical gradients can vary significantly within the footprint of a single satellite pixel. Beyond these radiometric uncertainties, thanks to their global coverage, satellite observations encompasses a wide range of optical water types and trophic states, exhibiting larger variability in both the environmental forcing and the phytoplankton community structure. The in situ BiOMaP dataset, in fact, is spatially concentrated in the Mediterranean Sea, the Black Sea, and adjacent Atlantic and Channel waters. While these waters do exhibit internal optical variability, they nonetheless occupy a relatively narrow and coherent region of the global bio-optical parameter space. The satellite matchup dataset used for fine-tuning, by contrast, spans the full range of oceanic provinces, from oligotrophic subtropical gyres characterized by extremely low particulate load, negligible CDOM, and absorption dominated by pure water, to coastal upwelling systems with high concentrations of both phytoplankton and non-algal particles, as well as polar waters, characterized by intense phytoplankton blooms accompanied by high concentrations of both non-algal particles and CDOM. These environments differ from the BiOMaP training domain not only in the magnitude of their optical constituents but also in the underlying covariance structure among the inherent optical properties (IOPs), posing a significant challenge for the CNN model. Finally, pigment concentrations included in the matchup dataset originate from multiple sources with non-uniform pigment suites and reporting conventions (e.g. naming or metadata), which further complicates the model training. Given these differences, the direct application of models trained on the regional in situ dataset to the global satellite one is not expected to yield optimal performance. Nevertheless, the previously trained models encode physically meaningful and biologically relevant correlations between spectral reflectance and pigment composition. To leverage this prior knowledge, a transfer-learning approach was adopted based on fine tuning. The weights of the best-performing CNN obtained in the previous phase were used as initialization, and the model was re-trained on the satellite-based dataset. Specifically, 100 epochs of learning were applied to the whole architecture, unfreezing all weights and using a reduced learning rate of 1 0 − 4 to preserve the previously learned representations while allowing adaptation to the new data distribution. To quantify the benefit of this strategy, the fine-tuned model (CNN F T) was compared to an identical architecture trained from scratch on the global satellite-based dataset, denoted as CNN SAT. In this latter configuration, all weights were randomly initialized and the full network was trained without constraints. This comparison allows to explicitly assess the added value of prior regional training and to evaluate the extent to which previously learned spectral–biogeochemical relationships improve convergence, stability, and predictive performance in a highly heterogeneous global context. Again, stratified splitting strategy into test/training set was applied, retaining 15% of the total sample as a test set and both noise injection and data augmentation techniques were considered for the training set. Specifically, standard deviations of 0.17 (17%) and 0.15 (15%) were applied for pigment concentration and satellite-based remote-sensing reflectance measurements, respectively, this latter supported by Zibordi et al. (2022)(corresponding to the average value of the satellite uncertainties excluding the longest center-wavelengths); and an augmentation factor of 100 were employed. Figure 8 compares the fine-tuned model (CNN F T) and the model trained from scratch (CNN SAT) using the same joint metric visualization applied in Figure 7: each symbol (bubble and square) corresponds to a pigment, the horizontal ad vertical axes indicate R 2 and M P E, while the symbol color represents the M A P E. Additionally, for each of the pigments, arrows connect the results obtained applying CNN SAT to those obtained applying CNN F T.

As expected, a general decrease of the model performance is noted with respect to those obtained on the in situ dataset. However, the fine-tuned model (CNN F T) consistently outperforms that trained from scratch (CNN SAT) across all considered pigments and all metrics, indicating that the prior knowledge acquired during training on the regional in situ dataset provides a measurable and systematic advantage when adapting the model to a global-scale satellite based application. Despite the increased radiometric uncertainty of the satellite observations and the possible spatial mismatch between satellite and HPLC measurements, the observed performance gains suggest that the convolutional layers of CNN F T have learned robust and transferable spectral features that are largely independent of the geographical location and representative of a wide range of optical water types. These features likely correspond to fundamental optical signatures associated with absorption properties of pigments included in the training dataset and consistently found across diverse oceanic regimes characterized by differing environmental conditions, phytoplankton assemblages, and observational uncertainties. The fine tuning strategy enables efficient adaptation to the new data distribution while mitigating the risk of overfitting to noise introduced by satellite-derived radiometry and heterogeneous in situ measurements. Improvements are observed not only for pigments that are optically dominant and well represented in the training data, but also for accessory pigments that typically exhibit weaker or more indirect spectral signatures. This result is particularly noteworthy given the increased uncertainty and missing-data patterns associated with the global HPLC compilation. It indicates that transfer learning enhances model robustness in data-scarce or noisy conditions, where training from scratch struggles to converge toward physically meaningful solutions. A complete view of the pigment-by-pigment metrics is given in Table 5.

TABLE 5. Performance metrics comparison on a pigment-by-pigment basis of CNN models trained based only on the satellite dataset (CNN SAT) and fine-tuned (CNN F T). Best values per metric are bolded.

The results indicate an average decrease, with respect to the in situ dataset based model, of approximately 0.18 in terms of R 2 after the fine-tuning (CNN F T vs. CNN), and approximately 0.29 applying the model trained from scratch (CNN SAT vs. CNN). Concerning M A P E and M P E, these generally increase by approximately 18% and 7% applying CNN F T, and by approximately 38% and 26% applying CNN SAT, respectively. Noteworthy is the retrieval of the concentration of allo, for which the fine-tuned satellite-based model slightly improves R 2 and substantially reduces relative errors with respect to the in situ CNN. This improvement is likely related to the average higher concentrations of allo represented in the satellite-based dataset with respect to the in situ dataset. In contrast, hex is among the pigments most negatively affected by the extension to satellite observations, showing a decrease of about 0.31 in R 2 and increases of approximately 21% and 23% in M P E and M A P E when applying CNN F T. In this case, however, no clear relationship can be identified between its representation in the input dataset and that in the satellite-based dataset.

3.3 Effect of spectral dimensionality

To further extend the applicability of the proposed model and to investigate the effect of reducing the input feature space, a reduced version of the CNN model, denoted as CNN *, was trained and evaluated. Selected wavelengths were 412, 442, 490, 555 and 670 nm, permitting the application to a number of legacy Ocean Colour missions as well as both the Copernicus-GlobColour and the European Space Agency’s Ocean-Colour Climate Change Iniciative (OC-CCI) datasets, thus supporting climatological studies. Pigment-by-pigment results are shown in Table 6.

TABLE 6. Performance metrics comparison on a pigment-by-pigment basis of CNN models trained with 11 (CNN) vs. 5 (CNN *) center-wavelengths. Best values per metric are bolded.

Reducing the spectral dimensionality leads to a limited but structured degradation in performance, confirming that a large fraction of the predictive skill is retained even under a strongly constrained spectral configuration. Overall, CNN * preserves the main performance hierarchy observed for the full CNN model, with changes in R 2, M S E and M A E remaining moderate for most pigments. In terms of explained variance ( R 2 ), differences between CNN * and CNN are generally small for the major pigments, with variations below or close to 0.02 for c h l i d e, c h l a, c h l b, c h l c 12, f u c o, d i a d, a l l o, z e a and c a r o. The largest degradation is observed for but, decreasing from 0.519 to 0.406, followed by hex and diato, both showing decreases of approximately 0.035. Conversely, peri shows a slight improvement when using the reduced spectral configuration, increasing from 0.839 to 0.863. Regarding the remaining metrics, similar patterns are obtained: CNN * generally exhibits slightly worse M S E, M A E and relative errors, although it improves selected metrics for some pigments, such as peri, allo, diad, and zea. Overall, the reduction to five center-wavelengths does not produce a systematic collapse in performance, but it does lead to more noticeable losses for pigments with weaker or more indirect spectral signatures, particularly but. The performance of CNN * is further summarized in Figure 9, using the same bubble representation adopted for the 11-center-wavelength configuration.

Complementing the quantitative performance comparison, a SHAP [SHapley Additive exPlanations, Štrumbelj and Kononenko (2013)] analysis was performed on both CNN and CNN * to identify the spectral bands driving pigment-specific predictions and to characterize how feature importance redistributes when the input dimensionality is reduced. SHAP analysis decomposes model predictions into individual feature contributions and is widely used in the remote sensing community to facilitate the interpretation of ML model outputs (e.g., Temenos et al., 2023; White et al., 2024; Marchese et al., 2026). Results for both configurations are shown in Figures 10a,b, respectively, where pigments and R r s ( λ ) are shown on the x- and y-axis, respectively, and the color indicates the mean absolute SHAP values expressed in percentage values relative to each pigment.

In both configurations, R r s ( 560 ) shows the highest relative contribution across most pigments. This is consistent with the role of R r s ( 560 ) as a well-established hinge point in c h l a retrieval algorithms, where green reflectance is particularly sensitive to variations in total pigment load. In most of the cases, the second and third highest contributors are R r s ( 442 ) and R r s ( 490 ). Both spectral regions coincide with the absorption features of several of the target pigments and are extensively used in existing c h l a and accessory pigment retrieval algorithms (Jeffrey, 2005; Roy et al., 2011; Clementson and Wojtasiewicz, 2019). In the 11-center-wavelength model, R r s ( 620 ) and R r s ( 681 ) emerge as additional significant contributors, particularly for c h l b, z e a, d i a t o, p e r i, b u t and h e x. These bands are absent from the 5-center-wavelength configuration, where their combined informational role is partially absorbed by R r s ( 673 ).

Similarly to the CNN model, CNN * was also fine-tuned based on a matchup dataset relying on the MAREDAT, NOMAD, and TARA datasets, and the Copernicus-Globcolour product consisting, in this case, of 7,287 observations. As for the 11-center-wavelength dataset, statistics highly differ from those of in situ dataset. Specifically, pigment concentrations show average values and standard deviations exceeding three and five times those included in situ dataset, respectively. Concerning remote-sensing reflectances, average values are on the order of those in the in situ dataset, however much higher variability is observed in this latter dataset, likely due to the marked presence of optically complex waters also reflected in the strong bimodal distributions obtained at 412, 442, and 490 nm (Supplementary Figure S7). Results obtained with both the fine-tuned version of CNN * (CNN F T *) and the CNN * re-trained from scratch (CNN SAT *) are summarized in Table 7.

TABLE 7. Performance metrics comparison on a pigment-by-pigment basis of CNN models trained based only on the satellite dataset (CNN SAT *) and fine-tuned (CNN F T *) with 5-center-wavelength. Best values per metric are bolded.

Analogously to the 11-center-wavelengths model, a general performance degradation is obtained with respect to the in situ dataset based model (CNN *) applying both CNN F T * and CNN SAT *. In particular, an average decrease of approximately 0.17 in R 2 is observed for CNN F T * and it reaches approximately 0.29 for CNN SAT *. Relative error metrics increase markedly, with average M A P E and M P E values reaching approximately 141% and 96% for CNN F T *, and approximately 184% and 137% for CNN SAT *, respectively. Nevertheless, the advantage of fine-tuning over training from scratch remains evident. This is especially clear in terms of R 2, where CNN F T * consistently achieves higher values across all 13 pigments, with significant improvements such as 0.61 versus 0.47 for chlide, 0.80 versus 0.73 for chla and 0.62 versus 0.50 for peri. The absolute error metrics further confirm the benefit of fine-tuning: CNN F T * systematically reduces M A E relative to CNN SAT * for all pigments, with particularly pronounced improvements for chla, fuco, chlc12, diad, and peri. Finally, the relative error metrics M A P E and M P E exhibit large values, particularly for pigments with very low concentrations such as chlide, zea, caro, and peri. Results are further illustrated in Figure 11, which compares CNN F T * and CNN SAT * using the same joint visualization of R 2, M P E, and M A P E.

4 Conclusion

This study evaluated the capability of simultaneously retrieving 13 diagnostic phytoplankton pigment concentrations from the multi-spectral band configuration of the OLCI sensor onboard Copernicus Sentinel-3. Five machine learning architectures, Random Forest, XGBoost, DNN, CNN, and BiLSTM, were systematically benchmarked under identical data and pre-processing conditions on a radiometrically rigorous in situ dataset spanning optically diverse European waters, constituting the first comprehensive ML architecture comparison specifically designed around the OLCI spectral configuration. The CNN architecture achieved the best overall performance, with coefficients of determination of 0.70–0.93 and mean relative errors of 25%–55% across the 13 target pigments. For the five best-retrieved pigments (chla, chlc12, fuco, caro, and diad) errors approach the approximate 20% inter-laboratory HPLC measurement absolute percent differences documented in Canuti (2023), indicating that the CNN results are near the physical limit of what can be expected from multi-spectral data for these compounds. The consistent outperformance of CNN over BiLSTM suggests that pigment retrieval from multi-spectral R r s is governed by local spectral features confined to a small number of consecutive bands, rather than long-range spectral dependencies, offering a recommendation for future architecture choices for ocean color retrieval problems. A novel two-phase Deep Learning training strategy was introduced and validated. By first optimizing pigment-specific modules independently before assembling them into a joint multi-output model with a shared dense head, the approach prevents the optimization dynamics of abundant, spectrally prominent pigments from dominating and degrading the retrieval of rarer accessory compounds. The extension to satellite observations was achieved through a transfer learning approach in which the CNN pre-trained on the high-quality regional BiOMaP in situ dataset was fine-tuned on a global matchup dataset combining MAREDAT, NOMAD, and TARA pigment measurements with Copernicus-GlobColour satellite Rrs products. A controlled comparison between the fine-tuned model (CNN F T) and an identical architecture trained from scratch on the satellite data alone (CNN SAT) provides a quantitative evidence that pre-training on radiometrically accurate regional in situ data confers a systematic and measurable advantage for global satellite pigment retrieval: R 2 decreased by approximately 0.2 for CNN F T versus 0.35 for CNN SAT relative to the in situ baseline, with corresponding improvements in MAPE between 20% and 40%. Specifically, CNN F T achieves average M A P E and M P E generally lower than 70% and 35%, compared to 100% and 50% for CNN SAT. Three pigment groups emerged with consistently distinct retrieval difficulty. Chla, fuco, chlc12, caro, diad, and peri are retrieved with R 2 larger than 0.85 and errors approaching the HPLC inter-comparison baseline, indicating very high maturity level. Chlide, chlb, hex, diato, and allo occupy an intermediate tier where DL techniques confers a meaningful advantage over classical ML approaches. Zea and but remain challenging regardless of model complexity, most likely due to low concentrations, poor representation in the training dataset, and high spectral overlap with dominant pigments, a behavior consistent with previously published studies (e.g., Pan et al., 2010; Bracher et al., 2015; Moisan et al., 2017; El Hourany et al., 2019; Stock and Subramaniam, 2020) and suggesting inherent limitation of multi-spectral retrieval that hyperspectral sensors such as PACE OCI may partially resolve. In addition to the OLCI-targeted model, a five-wavelength reduced version (CNN *) was developed and validated for compatibility with the Copernicus-GlobColour and ESA Ocean Colour Climate Change Initiative (OC-CCI) merged multi-sensor archives, enabling retrospective application of the multi-pigment retrieval framework to a continuous satellite record extending back to 1997. The 5-center-wavelength model exhibits similar in situ performance to the full 11-center-wavelength one for most pigments. For both models, the SHAP analysis indicates R r s ( 560 ), R r s ( 442 ), and R r s ( 490 ) as the mayor provider of informative content for most of the pigments, with additional significant contributions from R r s ( 620 ) and R r s ( 683 ) (or alternatively R r s ( 670 ) for the 5-center-wavelength model) for c h l b, z e a, d i a t o, p e r i, b u t and h e x. Similarly, the transfer learning advantage over training from scratch is preserved after fine-tuning.

It is recognized that for both the 11- and 5- center-wavelength models, the geographical concentration of the BiOMaP in situ training dataset in European waters may represent a limiting factor for global transferability; expanding coincident radiometric and HPLC measurements to underrepresented oceanic regions, including the tropical Pacific, Southern Ocean, and Atlantic and Pacific upwelling systems, would reduce the domain gap between regional pre-training and global satellite application and improve the robustness of transfer learning. Despite this limitation, this study demonstrates that CNN, combined with a two-phase training strategy, and a systematic transfer learning protocol, constitute a robust and operationally viable framework for multi-pigment retrieval from multi-spectral OLCI ocean color data. High-quality in situ radiometric measurements, free from the uncertainties introduced by atmospheric correction, adjacency effects, and sub-pixel heterogeneity, allows the model to identify spectral-biogeochemical relations that satellite-derived R r s alone cannot provide with the same accuracy, particularly for accessory pigments whose optical signatures are less evident and that are generally present with lower concentrations. Developed models are directly applicable to ongoing Copernicus Sentinel-3 mission’s and to the long-term Copernicus-Globcolor dataset. The applicability to the ESA OC-CCI dataset was not explicitly assessed, however no compatibility issues are expected. Additionally, the application to current and future hyperspectral missions including NASA PACE and ESA CHIME is envisaged. A large body of historical in situ HPLC measurements exists that was collected alongside multi-spectral radiometric observations spectrally incompatible with hyperspectral sensors, and would therefore be excluded from direct hyperspectral model training. The present results suggest that such datasets need not be discarded: provided their radiometric quality is sufficient, they can serve as pre-training foundations for hyperspectral CNN models, with the learned spectral-biogeochemical representations subsequently fine-tuned on the smaller but spectrally complete hyperspectral matchup datasets that PACE and CHIME will progressively generate.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. Generative AI tools were used solely to improve grammar, spelling, and overall readability, without altering the substantive content or interpretation of the text.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Sachidananda Mishra, National Oceanic and Atmospheric Administration (NOAA), United States

Reviewed by: Abdul Azeez Saleem, King Fahd University of Petroleum & Minerals, Saudi Arabia

Dilipkumar Jayaraj, UMR8187 Laboratoire d’océanologie et de géosciences (LOG), France