Hierarchical spectral ensemble with physics-informed augmentation for global water quality retrieval from hyperspectral remote sensing

1 Introduction

Inland and coastal water bodies are subject to increasing anthropogenic pressure. Eutrophication driven by nutrient loading, increased sediment flux from agricultural catchments and rising dissolved organic carbon inputs from terrestrial runoff all alter the optical and biogeochemical properties of water bodies faster than traditional monitoring networks can document (Diaz and Rosenberg, 2008; Chorus and Welker, 2021). The ecological consequences are large: decomposition of algal blooms causes hypoxic dead zones that now affect more than 400 coastal systems worldwide (Diaz and Rosenberg, 2008), while harmful cyanobacterial blooms pose direct risks to drinking water supply and aquatic biodiversity (Chorus and Welker, 2021). Understanding these changes at the global scale requires measurement approaches that are both spatially continuous and temporally frequent, requirements that cannot be met by in situ field campaigns given the approximately 304 million lakes and 37 million kilometres of rivers distributed across all continents (Verpoorter et al., 2014).

Satellite remote sensing offers the only operationally scalable pathway to global-scale water quality monitoring (Palmer et al., 2015; Sagan et al., 2020). The primary observable in passive aquatic remote sensing is the above-water remote sensing reflectance R _*rs* (λ), defined as the ratio of water-leaving radiance to downwelling irradiance at wavelength λ. The spectral shape of R _*rs* (λ) encodes the integrated optical contributions of all dissolved and particulate constituents in the water column through the radiative transfer relationship where a (λ) is the total absorption coefficient, b _*b* (λ) is the total backscattering coefficient, and f/Q accounts for the angular distribution of the underwater light field (Mobley et al., 2015; Odermatt et al., 2012). Because both a (λ) and b _*b* (λ) are additive functions of each optically active constituent, recovering concentrations from measured R _*rs* (λ) is fundamentally ill-posed when multiple constituents co-vary and their spectral signatures partially overlap.

Four ecologically relevant water quality parameters are targeted in this study. Chlorophyll-a (Chl-a, µg/L), the dominant phytoplankton pigment, is the most widely used indicator of eutrophication status and the risk of cyanobacterial blooms (Chorus and Welker, 2021), and most optical retrieval approaches are based on the absorption properties of Chl-a at 443 and 675 nm and the fluorescence emission at ∼683 nm. Total suspended solids (TSS, mg/L) control underwater light attenuation, sediment transport dynamics and benthic habitat quality (Wass et al., 1997); retrieval typically takes advantage of increased near-infrared (NIR) backscattering which increases with particle concentration (Nechad et al., 2010).

The absorption of coloured dissolved organic matter at 440 nm (aCDOM(440), m−1) is a measure of the contribution of terrestrially derived humic and fulvic acids to light absorption and has implications for dissolved organic carbon flux, carbon cycling and drinking water treatment (Brezonik et al., 2015). The Secchi depth (Z_sd, m) is a classical metric of water transparency and integrates the optical effects of all absorbing and scattering constituents and is related to the diffuse attenuation coefficient K_d through well-established mechanistic models (Lee et al., 2015; Lehmann et al., 2023).

Typical R_rs (λ) spectra of the five dominant water types in the GLORIA 2022 dataset are shown in Figure 1. The strong spectral differences between oligotrophic, mesotrophic, eutrophic, turbid and CDOM-rich regimes result from the bio-optical processes that are operating. Eutrophic spectra have a well-resolved fluorescence emission peak at about 683 nm and a phytoplankton red-edge at about 709 nm. Turbid waters have NIR backscattering enhancement beyond 750 nm. CDOM-rich waters have steeply decaying reflectance at short wavelengths. Crucially, the nearly identical spectral shapes for the low-concentration water types in the linear panel (Figure 1, left) demonstrate the basic inversion challenge in the oligotrophic regime, where the optical signals of constituents approach radiometric noise floors.

Simultaneous global retrieval of these four parameters from R _*rs* spectra is constrained by three structural factors that limit the performance of all existing approaches. First, target concentration ranges span four orders of magnitude (Chl-a: <0.05 to >500 μg/L; TSS: <0.1 to >2000 mg/L), producing highly right-skewed training distributions that cause standard regression estimators to optimise for the high-concentration tail while generating near-random predictions at low concentrations, even when reporting competitive aggregate metrics (Pahlevan et al., 2020). Second, the optical signatures of Chl-a, TSS, and a _CDOM(440) are not orthogonal: each constituent contributes to both absorption and scattering across overlapping wavelength regions, so the spectral contribution of any single constituent cannot be unambiguously isolated from measured R _*rs* (λ) (Odermatt et al., 2012). Third, the relationship between remote sensing reflectance Rrs(λ) and water quality parameters varies significantly across different ecosystems and depends upon the geographic regions. Variations in particle composition, phytoplankton communities, and dissolved organic matter sources make it hard for models trained on one regional dataset to perform accurately in unfamiliar or out-of-distribution environments (Pahlevan et al., 2020).

Early global water quality retrieval mainly depends on semi-analytical bio-optical models that calculates water properties by fitting absorption and backscattering coefficients to measured Rrs(λ) (Mobley et al., 2015; Odermatt et al., 2012). Although, these physics-based methods are easily interpretable and do not incline on the large labelled datasets, they are highly sensitive to assumptions regarding inherent optical properties, which changes considerably across different ecosystems, leading to declined performance in optically complex waters where multiple constituents changes simultaneously (Odermatt et al., 2012). The inventions of large and globally diverse hyperspectral datasets such as GLORIA (Lehmann et al., 2023) has introduced another promising direction: training empirical or hybrid models directly using co-located reflectance and water quality measurements. However, retrieval accuracy is still fundamentally limited by sensor spectral resolution, since multispectral satellite sensors with fewer than 10 bands cannot capture fine spectral characteristics, such as the Chl-a fluorescence peak near 683 nm or the NIR backscattering slope, that are needed for distinguishing water types under complex optical conditions (Palmer et al., 2015; Sagan et al., 2020), whereas airborne and in situ hyperspectral instruments with nearly 1 nm spectral resolution provide much richer spectral detail, though often with lower spatial coverage.

In recent years, machine learning techniques have shown clear improvements over traditional semi-analytical inversion models for getting all four water quality parameters (Sagan et al., 2020; Tian et al., 2023; Werther et al., 2022 as also highlighted by the paper Ma et al., 2023). Global-scale water quality analyser has improved considerably with the availability of multi-site and cross-ecosystem datasets, although accurately representing the entire concentration spectrum, from oligotrophic to hypereutrophic conditions, still remains difficult because most training datasets are dominated by mesotrophic and eutrophic samples (Pahlevan et al., 2020; Cao et al., 2020). Methods such as random forests (Soomets et al., 2020), gradient boosting machines (Tian et al., 2023; Werther et al., 2022), and mixture density networks (Pahlevan et al., 2020) have been successfully applied to inland and coastal water quality retrieval, showcasing that ensemble and probabilistic approaches are more capable of modelling the complex non-linear relationships between Rrs and constituent concentrations that conventional semi-analytical models often fail to capture. Uncertainty quantification has also become an very important aspect of retrieval evaluation, where probabilistic approaches such as MDN (Pahlevan et al., 2020) and Bayesian ensemble methods (Werther et al., 2022) have shown that estimating prediction uncertainty for individual samples provides more useful information than depending only on point estimates, especially under optically complex or data-sparse conditions. Importantly, uncertainty propagation throughout the complete retrieval chain all the way from sensor noise and atmospheric correction residuals to the regression model itself, has been recognised as significant for developing reliable operational products (Sagan et al., 2020; Werther et al., 2022), yet many existing studies still focus mainly on point-estimate accuracy without systematically analysing how uncertainty in the input data influences output uncertainty across different concentration ranges. Geographic transferability continues to remain a major challenge as it still being unresolved (Palmer et al., 2015; Sagan et al., 2020; Cao et al., 2020), since models that achieve strong performance within one ecosystem or region often experience a significant reduction in accuracy when applied to geographically independent water bodies, mainly because of changes in particle composition, phytoplankton communities, and dissolved organic matter characteristics (Odermatt et al., 2012). Spectral data augmentation has been explored as a possible solution for addressing class imbalance in skewed aquatic datasets (Pahlevan et al., 2020; Werther et al., 2022), although, rigorous leakage-free evaluation of augmentation effectiveness is still rarely reported in the literature. In addition to that, the geographic coverage, dataset composition, and validation strategies vary considerably across studies, making direct performance comparison inherently difficult (Palmer et al., 2015; Sagan et al., 2020).

Geographic transferability is a main challenge that differs from simple statistical overfitting. Even models trained using globally diverse datasets can perform poorly when applied to geographically independent water bodies, because regional variations in particle mineralogy, phytoplankton community composition, and dissolved organic matter characteristics can significantly change the relationship between Rrs and constituent concentrations in ways that a globally averaged model may not fully capture (Odermatt et al., 2012; Pahlevan et al., 2020). Werther et al. (Werther et al., 2022) demonstrated this issue in European lakes, where validation across different lakes produced considerably lower performance than validation within the same lake, mainly due to differences in optical water types between basins. Similarly, Sagan et al. (Sagan et al., 2020) identified geographic domain shift as one of the major problems preventing the operational deployment of machine learning-based retrieval models, signifying that successful training on data from one continent does not necessarily ensure good performance on another. With all of these concerns anyway, most published studies still rely on random train–test splitting approaches, which preserve geographic similarity between training and testing samples and therefore tend to overestimate real-world deployment accuracy (Palmer et al., 2015; Cao et al., 2020). More rigorous geographic hold-out evaluation strategies, where entire regions are excluded during training, are still relatively uncommon and remain absent from many comparative benchmarking studies.

Spectral data augmentation has been explored as a strategy to address class imbalance in skewed aquatic training distributions. Pahlevan et al. (2020) applied oversampling to under-represented oligotrophic bins in the MDN training pipeline, reporting marginal improvements in low-concentration retrieval. Werther et al. (2022) used bootstrap resampling and sample weighting to amplify the contribution of rare observations during ensemble training. However, most augmentation strategies reported in the aquatic remote sensing literature have been implemented prior to the train/test split, introducing the risk of information leakage: if synthetic spectra generated from training-set neighbours are present in the test partition, aggregate performance metrics may be inflated relative to true out-of-sample accuracy (Pahlevan et al., 2020; Werther et al., 2022). Rigorous evaluation of augmentation benefit under leakage-free conditions—where synthetic data are generated exclusively from within the training partition and withheld entirely from test evaluation—has not been systematically demonstrated for hyperspectral water quality retrieval.

Augmentation-based strategies to address oligotrophic data imbalance have received limited attention under rigorously leakage-free conditions, and per-concentration-stratum reporting remains non-standard despite the known masking effect of aggregate metrics on low-regime failure (Pahlevan et al., 2020; Werther et al., 2022). These gaps motivate the present study. Three systematic limitations additionally reduce the interpretability of performance metrics reported in the existing literature. Studies with R2 > 0.90 are mostly single-site assessments, where training and test data are from the same water body (Cao et al., 2020), and thus capitalize on uncontrolled geographic correlation between partitions, not a transferability to unseen environments. Beyond this, the near-universal practice of reporting only aggregate metrics is problematic for skewed distributions: because the squared-error contribution of high-concentration samples dominates the aggregate R ² computation, published values are disproportionately representative of eutrophic performance and can conceal systematic failure in the oligotrophic regime (Soomets et al., 2020; Mishra and Mishra, 2012). Another methodological issue is leakage of information from pre-processing, augmentation or feature-fitting operations performed prior to train/test partitioning, which may lead to overly optimistic test-set estimates that do not generalise to the true out-of-sample case.

In this paper, we address these limitations with four technical contributions: (1) the Hierarchical Spectral Ensemble, which fuses five heterogeneous base learners through a diversity-constrained NNLS meta-learner that outperforms individual models and simple baselines across all strata; (2) ISDA, an interpolation-based spectral augmentation approach with ecologically motivated bin boundaries applied only to training data after the stratified split, to avoid information leakage; (3) a 291-feature bio-optical engineering suite with all preprocessing, PCA fitting, and scaling operations performed only on training data; and (4) per-quartile stratified performance reporting that measures retrieval skill within each concentration stratum separately, exposing low-regime failures rather than masking them within an aggregate metric.

The specific objective of this study is to develop, validate and critically assess a globally trained machine learning framework for simultaneous retrieval of four ecologically critical water quality parameters from hyperspectral in situ Rrs spectra, with explicit attention to data leakage prevention, geographic generalisation, baseline comparison, feature explainability, and transparent reporting of performance limitations across the full concentration range. The term “globally trained” is used here deliberately because, while the framework is trained on a globally diverse dataset contributed by six continents, the geographic hold-out results illustrates that global training does not guarantee reliable global deployment, particularly for TSS and aCDOM(440) in spatially independent regions.

2 Dataset

2.1 GLORIA 2022

The Global Reflectance community dataset for Inland and Nearshore Areas (GLORIA) (Lehmann et al., 2023) is the largest publicly available global hyperspectral water quality database. The 2022 release includes 7,572 co-located R_rs spectra and in situ water quality measurements from over 1,300 unique water bodies across six continents. These water bodies represent ocean, coastal, estuarine, lacustrine, riverine, and reservoir environments. Each entry includes R _*rs* (λ) at 551 bands (350–900 nm, approximately 1 nm resolution), latitude, longitude, UTC timestamp, and laboratory measurements of Chl-a, TSS, a _CDOM(440), and Secchi depth.

The dynamic range column in Table 2 reports the ratio of the 99th-percentile to the first-percentile value of each measured parameter across the full GLORIA dataset prior to any outlier removal. For a _CDOM(440), values in GLORIA span from approximately 0.003 m⁻¹ in clear oceanic waters to approximately 30 m⁻¹ in heavily humic peat-draining rivers and dark-water tropical streams, yielding a raw dynamic range exceeding 10³×. This range is physically consistent with published field measurements in extreme peatland and blackwater river systems (Brezonik et al., 2015). After the 4.5×IQR outlier fence described in Section 2.2, the retained range reduces to the approximately 0.007–6.3 m⁻¹ span visible in Figure 2. Target availability and dynamic ranges are summarised in Table 2.

TABLE 1. Test-set performance: Condition A (no augmentation) vs. Condition B (post-split ISDA). All metrics evaluated on original, non-synthetic spectra only.

R ² and RMSE are computed in original physical units after back-transformation. 95% bootstrap confidence intervals from 1,000 paired resamples.

TABLE 2. Target variable availability and dynamic range in GLORIA 2022. Dynamic range is defined as the 99th/1st percentile ratio of the full pre-filtering dataset.

2.2 Preprocessing and splitting

Raw data undergo five preprocessing stages, with all statistical fitting operations performed exclusively on training data to prevent preprocessing leakage from validation or test partitions. Table 3 summarises the sample counts remaining after each filtering stage.

TABLE 3. Data filtering and preprocessing summary. All statistical fitting operations are determined from the training partition only.

Per-parameter sample counts differ because not all GLORIA, records contain measurements for all four target parameters (see Table 2). Test sets contain exclusively original, non-synthetic spectra in all experiments.

The 4.5×IQR outlier limitations is purposefully more lenient than the standard 1.5× threshold. The highly right-skewed distributions of aquatic parameters mean that a fence would wrongly remove valid hypereutrophic bloom events. Keeping concentration extremes is very important for a model trained on a globally diverse dataset which includes the entire trophic gradient. Savitzky-Golay smoothing (window = 11, polynomial order = 2) reduces the sensor noise while keeping the spectral shape and peak positions. Chl-a, TSS, and aCDOM(440) are transformed using log1p (·) to lessen skewness and stabilize regression variance. Zsd is kept in its original scale because its distribution is much less skewed. All R2 values reported in this paper are calculated in original physical units after back-transformation, applying exp (·) − 1 for the three log1p-transformed parameters unless noted otherwise. About 3.1% of the GLORIA 2022 records did not have a valid UTC timestamp. for these, a placeholder date of June 30 of the respective survey year was given to maintain temporal feature completeness without adding false precision.

3 Methods

3.1 Physics-informed feature engineering

From 551 smoothed reflectance bands, we constructed a 291-feature engineering suite across six groups. All feature engineering operations that involve parameter fitting are limited to training data. They are then only applied in transform-only mode to validation and test sets. Let R(λ) denote the smoothed reflectance.

• PCA spectral scores (80): Principal component analysis keeps 80 components, capturing more than 99% of cumulative variance. PCA offers a compact, noise-reduced representation of the full spectral shape. We fit the PCA model on training spectra only.
• Bio-optical spectral indices (52): These validated indices are based on established bio-optical theory and cover four parameter families. For Chl-a: NDCI (Mishra and Mishra, 2012), 3-band CI (Gitelson et al., 2008), 2-band ratio R (709)/R (665) (Gitelson et al., 2007), OC2/OC3/OC4 algorithms (O’Reilly et al., 1998), FLH (Gower et al., 2004), MCI (Matthews et al., 2012), CIre, and ChlGreen/ChlRed (Gitelson et al., 2003), along with UV/blue low-concentration indices. For TSS: NDSSI, Turb709 (Nechad et al., 2010), NIR/Red, NDTI, SPM, SABI, and low-turbidity indices. For aCDOM(440): SUV proxy, Bslope, UV-blue proxy ratios, and absorption slope indices (Brezonik et al., 2015). For Zsd: Kd proxy, Clarity index, and Kd (490) (Lee et al., 2015). Complete formulations of all spectral indices and engineered features are provided in Table 4.
• Derivative features (19): We use first and second Savitzky-Golay derivatives at seven key wavelengths (412, 443, 490, 560, 665, 709, 754 nm). Global derivative statistics can be able to record the spectral inflection points and rate-of-change information that band ratios alone will not be able to cover.
• Wavelet energy features (12): Five-level Haar decomposition gives multi-resolution spectral energy and mean coefficients per scale, which is able to capture both broad-scale spectral envelope shape and fine-scale features. Spatial features (26): Latitude and longitude are as raw values, squared terms, cross-products, and sinusoidal transforms. A K-means model (K = 25) fits training-set coordinates and assigns samples to 25 biogeographic optical water type zones through one-hot indicators.
• Spatio-temporal interactions (5):We use circular DOY encoding (sine/cosine over a 365-day period) interacted with latitude, longitude, and season to be able to capture the hemisphere-dependent phenological dynamics.

TABLE 4. Complete formulations of all 52 bio-optical spectral indices included in the HSE feature suite. R (λ) denotes R _*rs* (λ) in sr⁻¹. Indices are grouped by primary target parameter.

3.2 Interpolation-based spectral data augmentation (ISDA)

The chronic lack of oligotrophic samples in GLORIA, where Chl-a below 1 μg/L accumulates for about 8.6% of records, leads standard regressors to give out almost random predictions in the low-concentration range, even though they report competitive aggregate R2 values. ISDA tackles this by creating synthetic training spectra through convex interpolation between neighboring real spectra within ecologically defined concentration bins. The procedure is outlined in Algorithm 1.

Bin edges follow established ecological trophic thresholds (Table 5), ensuring interpolation occurs between spectrally similar water types rather than across distinct ecological regimes. The Beta(0.4, 0.4) prior is bimodal, steering synthetic spectra toward real endpoints to keep spectral realism intact. We apply ISDA only to the training set, after the stratified 70/15/15 train/validation/test split. Validation and test sets do not include synthetic samples in any experiment, ensuring test-set performance reflects model behavior on real spectra only. Table 5 summarizes the pre- and post-augmentation training set sizes.

TABLE 5. ISDA ecological bin boundaries and pre/post-augmentation training set sizes (Condition B, post-split augmentation).

Synthetic percentages are intentionally modest; ISDA, targets only the most severely under-represented oligotrophic bin. The five-times cap (C = 5.0) prevents synthetic samples from dominating any individual bin.

Algorithm 1 ISDA: Ecologically motivated spectral augmentation (training-set only).

3.3 Hyperparameter optimisation and ensemble learner selection

Five base learners were selected to maximise ensemble diversity while ensuring complementary bias-variance profiles. Starting from a candidate pool of seven tree-based algorithms (XGBoost (Chen and Guestrin, 2016), LightGBM (Ke et al., 2017), HistGBM, CatBoost (Prokhorenkova et al., 2018), ExtraTrees (Geurts et al., 2006), Random Forest, and Gradient Boosting), models were added greedily in order of their marginal contribution to ensemble R ² on the validation set. Random Forest and Gradient Boosting were excluded from the final ensemble because their marginal contribution was below 0.002 R ² units (mean across all four targets) while adding approximately 35% to total training time.

Optuna (Akiba et al., 2019) with TPE sampling tunes XGBoost and LightGBM (60 trials each) and ExtraTrees (30 trials) against the held-out validation set. The search covers tree depth (3–10), learning rate (5 × 10⁻³–0.15, log-uniform), regularisation (10⁻⁴–10, log-uniform), and subsampling fractions (0.5–1.0). The optimal tree count found during HPO defines a fixed upper bound n _est = 1,200 used in the OOF loop, where no early stopping is applied.

3.4 Out-of-fold stacking

Ten-fold stratified OOF stacking generates M ∈ R^*n*tr×5. Each model trains on nine folds and predicts on the held-out fold, iter-ated over all 10 splits. The five base learners are: (1) XGBoost (Chen and Guestrin, 2016): 1,200 trees, pseudo-Huber objective; (2) LightGBM (Ke et al., 2017): 1,200 trees, Huber (α = 0.9); (3) ExtraTrees (Geurts et al., 2006): 800 trees; (4) HistGBM: 800 iterations, 255 bins; (5) CatBoost (Prokhorenkova et al., 2018): 1,200 trees, ordered boosting, Huber loss. Sample weights inversely proportional to bin frequency are applied during OOF training to amplify the loss contribution of rare low-concentration samples.

3.5 Diversity-constrained NNLS meta-learner

Non-negative least squares (NNLS) regression combines the five sets of OOF predictions into the final ensemble. NNLS ensures non-negative combination weights without explicit weight clipping, providing convex-combination properties appropriate for ensemble blending. The regression intercept is retained for bias correction. In validation experiments, NNLS outperformed constrained Ridge regression; Ridge degraded Chl-a and TSS R ², consistent with the observation that Ridge can assign positive weight to base learners whose OOF predictions are negatively correlated with the target in a given fold, whereas NNLS constrains the search to non-negative contributions.

3.6 Baseline models

Four baseline models are trained on the same preprocessed features and train/test split to provide a systematic reference: Ridge regression, Decision Tree (maximum depth 5), Random Forest (100 trees, maximum depth 15), and XGBoost with default hyperparameters. No hyperparameter optimisation, stacking, or ISDA augmentation is applied to baselines.

The HSE pipeline architecture is illustrated in Figure 3.

4 Experimental results

4.1 Leakage-free ISDA and condition A/B comparison

To directly measure the impact of augmentation while ensuring that no information leakage occurs, two experimental conditions are considered. Condition A uses a standard stratified 70/15/15 data split without including any synthetic samples in the dataset. Condition B applies ISDA only to the training partition after the split, while the validation and test sets remain completely unchanged. Both experimental settings follow the same preprocessing pipeline, hyperparameter optimisation strategy, and stacking methodology, ensuring that any observed performance difference can be attributed entirely to the use of ISDA.

Post-split ISDA (Condition B) results in a modest improvement for Chl-a (+0.007 R2) and aCDOM(440) (+0.005 R2), while showing almost no change for Zsd and no meaningful improvement for TSS (−0.002 R2). The highly similar outcomes observed between Conditions A and B indicate two important findings: first, the leakage-free implementation does not artificially increase overall test-set performance when compared with the no-augmentation baseline, and second, ISDA contributes only a limited improvement in overall aggregate R2 performance. This is expected, given that the available oligotrophic pool in GLORIA is too small for interpolation to meaningfully bridge the spectral indistinguishability documented in Table 6. Figure 4 visualises the comparison, and Table 1 presents full numeric results.

TABLE 6. Base learner OOF R ² and HSE stacked result, evaluated in log_1*p*-transformed target space (Condition B). OOF R ² in log space differs from back-transformed test-set R ² in Table 4 due to distributional compression at high concentrations.

4.2 Test-set performance, per-quartile analysis, and ensemble composition

Figure 5 shows the predicted versus measured scatter for all four parameters on the held-out test set (Condition B). HSE achieves R ² ≥ 0.84 for a _CDOM(440) and Z _*sd*, and R ² = 0.822 and 0.704 for Chl-a and TSS respectively, with a mean R ² = 0.832 across all four parameters.

Although these overall results appear encouraging, they hide a more complex underlying trend. Table 6 shows the concentration-stratified Q1–Q4 breakdown for Condition B. For Chl-a, TSS, and aCDOM(440), the R2 values across Q1–Q3 are negative, indicating that HSE performs worse than a simple mean-based predictor within the low-to-moderate concentration ranges. As a result, the overall aggregate R2 values are driven almost entirely by the strong performance observed in the high-concentration Q4 range. In contrast, Zsd demonstrates positive R2 values from Q2 onward, likely due to its relatively smaller dynamic range and the more consistent spectral relationship between water transparency and diffuse attenuation. Taken together, these findings show that the strong aggregate R2 performance of the framework does not fully reflect its behaviour across all concentration levels.

Masks an almost complete failure in predictive performance across most of the concentration range for three of the four target parameters. Several observations from the per-quartile analysis require further interpretation. The strongly negative R2 value in Chl-a Q2 (−1.999) is especially significant because Q2 represents the 1.72–4.70 μg/L mesotrophic transition range, where optical signals from water constituents are present but not sufficiently dominant. This creates ambiguity in the spectral response, causing the ensemble predictions to deviate systematically from the measured values and indicating that spectral information alone is inadequate for reliable retrieval within this range. In the case of TSS, the Q4 R2 value of only +0.121, despite an overall aggregate R2 of 0.704, highlights how heavily the aggregate metric is influenced by the high-concentration tail; even at elevated concentrations, TSS retrieval performance remains below the reliability expected for many operational monitoring applications. For aCDOM(440), negative R2 values continue across Q1–Q3, reflecting the well-known spectral overlap between CDOM absorption and chlorophyll-a absorption within the blue wavelength region (Odermatt et al., 2012; Brezonik et al., 2015). This type of spectral ambiguity has been widely reported in both field-based bio-optical studies (Brezonik et al., 2015) and large multi-site machine learning evaluations (Sagan et al., 2020; Werther et al., 2022), where separating the CDOM signal from phytoplankton contributions becomes particularly difficult at low concentrations. Zsd is the sole parameter exhibiting positive R2 from Q2 upward; its physically motivated relationship to the diffuse attenuation coefficient (Lee et al., 2015) provides a more stable spectral signal that is less dependent on constituent-specific optical properties. These parameter-level differences in per-quartile behaviour point to distinct retrieval challenges that aggregate metrics entirely obscure, and they should inform the design of targeted future retrieval strategies for each parameter separately. The residual distributions for all four target variables are shown in Figure 6, demonstrating that the prediction errors remain centered around zero with no substantial systematic bias across the prediction range. The practical implication is that any operational use of HSE must account for which concentration regime the target water body falls in: the framework is reliable only for high-concentration eutrophic and turbid waters, and outputs in oligotrophic regimes should be treated as unreliable without additional uncertainty quantification.

Table 7 reports the OOF R ² of each base learner in log_1*p*-transformed target space. The stacked ensemble consistently improves over the best individual model by 0.001–0.012 R ² units, confirming that the NNLS meta-learner exploits non-redundant predictive information across the five base learners. This marginal but consistent improvement across all four parameters supports the ensemble design, even though the gains are modest relative to the best individual learner. Figure 7 presents feature importance scores for all four parameters.

TABLE 7. Base learner OOF R ² and HSE stacked result, evaluated in log_1*p*-transformed target space (Condition B). OOF R ² in log space differs from back-transformed test-set R ² in Table 1 due to distributional compression at high concentrations.

4.3 Baseline comparison and literature benchmarking

Table 8 quantifies the performance improvement attributable to the full HSE framework relative to the four baseline models. HSE produces consistent R ² improvements for Chl-a (+0.109 over Ridge) and a _CDOM(440) (+0.141 over Ridge). The Ridge regression result of R ² = −0.410 for TSS is itself informative: it confirms that the R _*rs*–TSS relationship is strongly non-linear, and that linear models are inadequate even as rough baselines for this parameter. Figure 8 illustrates the comparison.

TABLE 8. Baseline comparison: HSE ensemble vs. four simple baseline models. No hyperparameter optimisation, stacking, or ISDA is applied to baselines.

R ² values are in original physical units after back-transformation. ΔR ² denotes HSE minus baseline model.

Bold values indicate the best-performing result for each metric.

Table 9 benchmarks HSE against published algorithms. These comparisons are necessarily indicative rather than definitive, because direct cross-study comparison is confounded by differences in datasets, preprocessing procedures, sensor characteristics, and validation frameworks. Studies based on single lakes or single ecosystems (Cao et al., 2020) are not directly comparable with a global framework. Regarding the MDN Chl-a result (R ² = 0.37) attributed to Pahlevan et al. (2020): this value is derived from applying the MDN algorithm to the GLORIA 2022 test set used here, and is therefore a re-evaluation on a different dataset rather than a reproduction of published results.

TABLE 9. Comparison of HSE with published water quality retrieval methods.

Single-site or single-ecosystem evaluation.

Re-evaluated on GLORIA 2022 test set using back-transformed R2; original study reports different metrics.

5 Discussion

5.1 Global versus site-specific performance

The overall HSE performance obtained under random splitting (R2 = 0.822–0.962) is comparable to many published global retrieval algorithms tested on similar datasets. However, an important limitation must be considered: random-split evaluation still preserves geographic correlation between the training and test datasets, since GLORIA samples collected from geographically nearby water bodies are likely to be distributed across both partitions. The geographic validation results demonstrate substantial performance degradation on spatially independent regions, confirming that random-split aggregate metrics overestimate the geographic generalisation capability of the framework. These results should therefore not be interpreted as evidence of global applicability across all regions and concentration ranges, and we present the geographic hold-out figures in Table 9 alongside the random-split results to make this limitation explicit.

It is also worth noting that near-perfect accuracy on a globally heterogeneous dataset would itself be a warning sign of overfitting rather than genuine generalisation. Residual prediction variance is partly irreducible: it reflects in situ measurement uncertainty of approximately 10%–20% replicate variability for Chl-a (Lehmann et al., 2023), atmospheric correction residuals, and the intrinsic optical complexity of a globally diverse dataset. These sources of error are present regardless of model design.

5.2 TSS retrieval and spectral equifinality

The TSS R ² = 0.704 under random splitting, the deeply negative R ² values in Q1–Q3 (Table 6), and the geographic hold-out failure (Table 9) each reflect distinct aspects of the fundamental difficulty of passive optical TSS retrieval. At the core of this difficulty is what might be called spectral equifinality: mineral particles, biogenic detritus, and flocculated aggregates can produce nearly identical R _*rs* spectra despite mass concentrations differing by an order of magnitude (Nechad et al., 2010). This is not a modelling failure that more training data or better feature engineering can fully resolve, it represents a physical information ceiling that passive optical methods cannot overcome without additional constraints on particle composition or size distribution. The implication is that operational TSS retrieval from multispectral or even hyperspectral imagery will require either auxiliary information about particle type, or a fundamental shift toward active sensing approaches for high-confidence retrievals. This spectral equifinality problem is not unique to the present study; it has been documented in field-based bio-optical investigations (Odermatt et al., 2012; Nechad et al., 2010) and in multi-site machine learning evaluations (Sagan et al., 2020; Tian et al., 2023), where TSS consistently exhibits larger prediction uncertainty and weaker cross-site transferability than Chl-a or Zsd. Approaches incorporating prior information on particle size distribution or mineralogy have been proposed as a path forward (Odermatt et al., 2012), but such auxiliary data are rarely available at global scale.

5.3 Geographic generalisation and spatial feature dependence

The geographic validation experiments reveal a design tension in the current HSE framework. Including spatial features (latitude, longitude, K-means biogeographic zone indicators, and their interactions) improves random-split performance by enabling the ensemble to learn region-specific R _*rs*–concentration relationships. However, these same features degrade geographic generalisation because they encode location-specific optical information that does not transfer to spatially independent test regions. For Z _*sd*, geographic generalisation is comparatively robust, likely because water transparency has a more universal spectral relationship driven by total diffuse attenuation that is less sensitive to region-specific optical water types. Practitioners intending to apply HSE to novel geographic areas should treat the spatial features with caution and may benefit from retraining with locally available data. The geographic transferability gap observed here—particularly the TSS hold-out R2 of −0.200 and aCDOM(440) hold-out R2 of 0.152, is consistent with the broader literature on geographic domain shift in aquatic machine learning models (Sagan et al., 2020; Pahlevan et al., 2020), where models trained on globally diverse data nonetheless fail in spatially independent regions due to differences in particle composition, dissolved organic matter sources, and regional optical water type distributions (Odermatt et al., 2012). This highlights that although global training coverage is important, it alone is not enough to ensure reliable geographic generalisation, and therefore the random-split aggregate R2 values reported in Table 1 should be interpreted as upper-bound estimates of framework performance rather than true indicators of deployable accuracy. The variation in transferability observed across different parameters, where Zsd generalises relatively well while TSS performs very poorly, reflects the physical difference between parameters that possess a more universal spectral relationship, such as water transparency linked to diffuse attenuation, and those that depend heavily on region-specific particle composition or dissolved matter characteristics. This parameter-specific variation in geographic transferability is itself an important finding and may help guide future decisions related to sensor selection, validation strategies, and the development of regionally stratified retrieval frameworks.

5.4 The oligotrophic regime

The per-quartile analysis in Table 6 identifies what we regard as the primary limitation of the current framework, and one that applies broadly to data-driven retrieval methods. Negative R ² values in Q1–Q3 for Chl-a, TSS, and a _CDOM(440) demonstrate that the aggregate R ² values substantially misrepresent retrieval performance: the framework achieves reliable predictive skill only in the high-concentration Q4 regime. This is not simply a consequence of insufficient training data in oligotrophic bins. In oligotrophic waters, constituent optical signals approach radiometric noise floors, reducing the signal-to-noise ratio of the spectral input below the discrimination threshold of any data-driven retrieval model (Pahlevan et al., 2020). Compounding this, in situ measurement uncertainty of approximately 10%–20% for low-concentration samples (Lehmann et al., 2023) establishes an irreducible lower bound on achievable R ² in these regimes, independently of model design.

The results suggest two specific priorities for future work. First, physical-based synthetic spectrum generation with radiative transfer simulations, e.g., Hydrolight (Mobley et al., 2015) could generate physically consistent oligotrophic spectra with spectral diversity that is not possible with interpolation; this approach has been demonstrated in adjacent coastal remote sensing contexts (Odermatt et al., 2012) and would directly address the spectral indistinguishability documented in Table 1. Second, the use of uncertainty-aware training objectives that account for heteroscedastic measurement error in ground-truth labels can potentially reduce the adverse effect of label noise on low-regime predictions, consistent with recent probabilistic retrieval approaches (Werther et al., 2022) that quantify per-sample prediction uncertainty alongside point estimates. More broadly, a complete uncertainty budget for operational water quality retrieval should characterise how uncertainty propagates from sensor noise and atmospheric correction residuals through to model predictions; this is an under-studied problem that frameworks such as MDN (Pahlevan et al., 2020) and Bayesian ensemble methods (Werther et al., 2022) have begun to address, but which requires explicit treatment across the full concentration range, particularly in oligotrophic regimes where signal-to-noise is lowest. The observed oligotrophic failure is consistent with findings reported for the MDN retrieval framework (Pahlevan et al., 2020), which also demonstrated diminished performance in low-concentration regimes on GLORIA, and with the broader machine learning water quality literature in which Q1 performance is rarely reported and aggregate metrics dominate (Sagan et al., 2020; Werther et al., 2022; Soomets et al., 2020). Until such or similar advances are made, the reporting of aggregate R2 values, along with per-quartile breakdowns, should be considered standard practice in the field.

5.4.1 Sustainability implications

The increasing rate of global eutrophication (Chorus and Welker, 2021), the expansion of hypoxia zones (Diaz and Rosenberg, 2008), and the intensification of harmful algal blooms underline the need for operational water quality monitoring at scales that cannot be reached with field campaigns. The HSE pipeline offers a transferable, open-source methodological baseline whose design principles—leakage-free augmentation, per-quartile reporting, and diversity-constrained ensemble learning—may inform the development of future satellite-derived water quality products from Sentinel-3 OLCI, Sentinel-2 MSI, MODIS, and the upcoming PACE OCI hyperspectral imager. It should be emphasised, however, that the present framework operates on in situ hyperspectral Rrs and has not been validated on satellite-derived imagery; the path to operational satellite deployment remains substantially open. However, transition from this research framework to operational deployment would require substantial additional steps: spectral resampling and retraining for specific satellite sensors, thorough validation on independent datasets with sensor-specific atmospheric correction, and demonstrated performance in oligotrophic and spatially novel environments where the current framework shows the most significant limitations. As the geographic hold-out and per-quartile results confirm, the framework in its current form should not be applied operationally without such adaptations, particularly for TSS and aCDOM(440) retrieval. The primary contribution of this work for sustainability monitoring is therefore methodological: it provides a validated, transparent, and reproducible framework that demonstrates both the current capabilities and the current limits of data-driven hyperspectral water quality retrieval, with direct relevance to the goals of UN SDG 6 and SDG 14. Further validation and sensor-specific adaptation remain prerequisites before operational deployment can be recommended.

6 Limitations

The following limitations are acknowledged to support transparent interpretation of the reported results.

Oligotrophic retrieval failure. The model produces R ² < 0 in Q1–Q3 for Chl-a, TSS, and a _CDOM(440), emphasizing predictions less accurate than a mean baseline within these strata. Aggregate R ² values are driven by high-concentration Q4 performance and are not representative of retrieval accuracy in oligotrophic or low-to-moderate concentration regimes. This is the most consequential limitation of the current framework.

Geographic generalisation. Geographic hold-out experiments demonstrate substantial performance degradation for TSS and a _CDOM(440) on spatially independent test regions. The spatial features included in the framework should be treated with caution in deployment to novel geographic areas.

Hyperspectral input requirement. All experiments use in situ hyperspectral R _*rs* at approximately 1 nm resolution. The model cannot be applied to multispectral satellite imagery without spectral resampling, retraining, and adaptation for atmospheric correction artefacts.

Instrument and measurement variability. GLORIA 2022 aggregates measurements from multiple field campaigns using different radiometers and laboratory protocols. Inter-instrument systematic differences have not been independently harmonised, contributing to the effective noise floor in both training and evaluation.

Random-split versus geographic-split performance gap. The reported aggregate R ² values are from random stratified splits that retain geographic correlation. Geographic hold-out results are substantially lower and represent a more conservative assessment of generalisation.

Augmentation scope. ISDA is an interpolation-based statistical augmentation strategy; the ecological motivation refers to bin boundary selection rather than to physically constrained spectrum synthesis. Augmentation effectiveness is limited by the small oligotrophic source spectral pool in GLORIA, and the marginal aggregate gains in Table 1 reflect this constraint.

7 Conclusion

This paper presented the Hierarchical Spectral Ensemble (HSE), a leakage-free machine learning framework for water quality parameters retrieval from hyperspectral in situ remote sensing reflectance. Applied to GLORIA 2022 (7,572 samples, 551 bands, six continents), HSE achieves R ² = 0.822, 0.704, 0.841, and 0.962 for Chl-a, TSS, a _CDOM(440), and Z _*sd* respectively under leakage-free random-split evaluation (mean R ² = 0.832; all values in original physical units after back-transformation).

Four methodological contributions distinguish this work: (1) diversity-constrained NNLS meta-learning over five hetero-geneous base learners; (2) post-split ISDA augmentation with ecologically motivated bin boundaries applied exclusively to training data; (3) a 291-feature bio-optical engineering suite with complete leakage-free preprocessing; and (4) per-quartile stratified performance reporting that makes concentration-regime failures visible rather than absorbed into aggregate metrics.

Systematic baseline comparisons confirm consistent R ² improvements over Ridge regression, Decision Tree, Random Forest, and default XGBoost, with the largest gains for Chl-a (+0.109 R ² over Ridge) and a _CDOM(440) (+0.141 R ² over Ridge).

The per-quartile analysis leads to three principal findings that qualify these results. First, aggregate R ² values conceal near-complete predictive failure in oligotrophic Q1–Q3 strata for Chl-a, TSS, and a _CDOM(440); reliable predictive skill is achieved only in high-concentration Q4 regimes. Second, geographic hold-out experiments reveal substantial generalisation failure for TSS and a _CDOM(440) on spatially independent regions, indicating that the reported aggregate R ² values overestimate deployable geographic transferability. Third, the leakage-free ISDA implementation produces only marginal aggregate improvements (ΔR ² ≤ 0.007), because the available oligotrophic spectral pool in GLORIA is insufficient to address the fundamental spectral indistinguishability in the low-concentration regime. Addressing these limitations, particularly through radiative transfer-based synthetic spectrum generation and uncertainty-aware training objectives, represents the most important direction for future work.

All code and trained model artefacts are publicly available to support reproducible water quality monitoring. This work contributes to the methodological foundation for scalable water quality assessment relevant to SDG 6 and SDG 14, with the explicit acknowledgement that current oligotrophic regime performance and geographic generalisation capability require substantial further development before operational deployment can be recommended.

8 Reproducibility statement

All experiments were implemented in Python 3.10. Hardware: 2 NVIDIA GPU (CUDA-enabled), 64 GB RAM, 16-core CPU. Key dependencies: XGBoost 1.7 (Chen and Guestrin, 2016), LightGBM 3.3 (Ke et al., 2017), CatBoost 1.2 (Prokhorenkova et al., 2018), scikit-learn 1.2, Optuna 3.0 (Akiba et al., 2019), NumPy 1.23. Random seed: 42, fixed in NumPy (np.random.seed (42)), scikit-learn (random state = 42), and all tree-based models via their respective random state parameters; stratified splits are deterministic given this seed. Wall-clock runtime: approximately 4–6 h per target parameter for the full pipeline including Optuna HPO (60 trials for XGBoost/LightGBM, 30 for ExtraTrees), OOF stacking, geographic validation, and zone cross-validation on the specified hardware. All code and trained model artefacts are publicly available.

Acknowledgements

The authors gratefully acknowledge Sri Mata Amritanandamayi Devi (Amma), Chancellor, Amrita Vishwa Vidyapeetham, for her inspiration and for providing financial support for the Article Processing Charges (APC) of this publication.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

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: Pedzisai Kowe, Midlands State University, Zimbabwe

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

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