## 3 Results
### Correlation and Regression Analyses  
Correlation analysis was applied to examine relationships between in situ Chl-\(a\) concentrations and the corresponding spectral indexes derived from atmospherically corrected matchup pixels [84; 85]. Figure 5 presents the resulting correlation coefficients for the retrieval indexes. A heatmap of only R\({}^{2}\) values from the correlation analysis is also provided in Figure 52 of the Supplementary Materials. These figures reveal that some indexes, such as I\({}_{27}\), exhibited an overall negative correlation, that is, an increase in the index value corresponds to a decrease in Chl-\(a\) concentration. Conversely, indexes such as for example I\({}_{26}\), exhibited a positive correlation for most scenarios, indicating a direct relationship between the index value and the Chl-\(a\) concentration. However, the majority of correlation metrics were in the range \(\pm 0.5\), indicating statistically insignificant correlation. Furthermore, the results indicated that, regardless of the data category, both Sentinel-2 and Landsat 8 constantly had superior correlations compared to the other satellites. Of note is that while correlation analysis shows the strength and direction of the relationship, it does not provide insights into the models' prediction accuracy or biases. To select the optimal scheme (combination of atmospheric correction processor and Chl-\(a\) retrieval index) for each scenario (data category and satellite), we primarily relied on the RF importance scores because they offer a more robust feature selection metric for regression models compared to traditional correlation coefficients. Therefore, for each scenario, we initially ranked schemes by RF importance score. The top-ranked scheme was then selected, unless the second-ranked scheme demonstrated significantly higher correlation coefficients (for both Pearson's \(r\) and Spearman's \(\rho\)). Out of 28 scenarios, the optimal schemes of 20 scenarios were directly determined by the highest RF importance score, while the remaining 8 had the second-highest RF score but ranked significantly better based on the correlation coefficients. Next, a linear regression was fitted between the scheme-derived and matchup Chl-\(a\) concentrations to generate a predictive equation for Chl-\(a\). The predicted versus matchup (measured) Chl-\(a\) concentrations are plotted in Figure 6. Details on the performance metrics as well as the regression equations, are presented in Table S3 of Supplementary Materials. The plots in Figure 6 are colour-coded by satellite (Landsat 5, 7, 8, and Sentinel-2, left to right) and identify the four in situ data categories defined in Section 2.2: all (no filtering), location (HH or WLO), seasonality (AW or SS), and TSI (OM or EH). The retrieval index and atmospheric correction processor as well as the values for the matchup count, regression slope, RMSE, and R\({}^{2}\) are also shown on the plots. Figure 4: Random Forest feature importance analysis with colour-coded atmospheric corrections processors. The \(x\)-axis denotes the retrieval index (feature), and the \(y\)-axis shows the importance score (unitless). For each scenario, the most significant scheme is marked with an asterisk. Figure 5: Evaluation of schemes across different scenarios based on correlation analysis. Marker colours denote different atmospheric corrections, shapes represent satellites, and sizes signify the number of matchups. Landsat 5 and 7 data typically exhibited higher RMSE, from 1.5 (WLO) to 14.1 \(\upmu\)g/L (EH). In contrast, Landsat 5's HH and EH categories had nearly zero slopes, indicating ineffective predictive power. The best performance for Landsat 5 occurred in the SS category,  
Figure 6: Plots comparing modeled vs. measured Chl-\(a\) concentrations across satellites (columns) and data categories (rows), demonstrating the regression models’ performance. with a slope and R\({}^{2}\) of 0.38, indicating overprediction and underprediction at low and high concentrations, respectively.