Remote Sensing and Machine Learning for Ammonium (NH_4⁺) Estimation in the Oligotrophic Coastal Aquaculture Waters of Pagasitikos Gulf (Eastern Mediterranean)

1. Introduction

Intensive fish farming takes place in floating sea cages which are placed near the surface of the water column, particularly in coastal waters. It is characterized as an open cultivation system, where water passes freely through the cages, and interacts with the marine environment by producing effluents that are released directly into water. The metabolic processes of cultivated fish and the uneaten feed tend to result in nutrient release [1]. Although marine cage aquaculture contributes only 1% to total nutrient enrichment in coastal water bodies compared with other anthropogenic activities [2], it remains a relevant source of nutrient loading particularly in semi-enclosed areas.

Ammonium (NH₄⁺) is the main ionized nitrogen (N) form produced by fish farming [3] and it is the most preferred N source for phytoplankton growth [4]. Therefore, systematic monitoring of NH₄⁺ in areas where fish farming takes place is essential, as increased NH₄⁺ flux may stimulate excessive phytoplankton growth, which, under certain environmental conditions, could affect the ecological status of the surrounding ecosystem. In aquaculture environments, the estimation of water quality parameters, including NH₄⁺, is particularly based on field measurements, sampling and laboratory analysis. Although these methods provide extremely accurate results, they are costly, time-consuming and they are only able to provide point estimations. Recently, remote sensing has been applied as an alternative approach for estimating water quality parameters (WQPs) in both coastal [5,6,7,8] and inland [9,10,11,12] water bodies, associated with aquaculture operations, providing cost-effective, large-scale datasets with high temporal frequency [13].

Unlike optically active constituents (OACs) of water, the estimation of NH₄⁺ using remote sensing remains challenging because NH₄⁺ lacks a direct optical signature and does not significantly influence the absorption and scattering properties of water. In general, the retrieval of nutrients from remotely sensed data has relied on empirical, machine learning (ML) and deep learning (DL) approaches [14], with model development primarily exploiting the indirect relationships between non-OACs and OACs, such as chl a, total suspended matter (TSM) and colored dissolved organic matter (CDOM) [15]. However, these relationships may vary across water bodies due to differences in spatial and temporal characteristics, as well as trophic conditions, which regulate OAC abundance and consequently affect optical properties and the retrieval of non-OACs [14].

An alternative approach involves the direct estimation of non-OACs from remotely sensed data without prior retrieval of OACs [16,17]. This strategy relies increasingly on ML and DL approaches capable of capturing complex non-linear relationships between spectral bands and non-OACs of water [14,17]. Comparative studies have shown that ML models may demonstrate better generalization performance than DL models under conditions of limited data availability, as they are generally less prone to overfitting when trained on small datasets [14]. Overfitting is one of the major challenges in ML because models that perform well on training datasets may fail to generalize to unseen data [18].

The majority of remote sensing-based studies on nutrient estimation has focused on eutrophic or hypertrophic water bodies influenced by agricultural, industrial, and/or urban wastewater runoff, while the effect of seasonal variability has received limited attention [14]. This study aims to investigate the potential of multispectral Sentinel-2 MSI imagery for estimating NH₄⁺ in an oligotrophic semi-enclosed cove where a marine cage aquaculture operates. Firstly, the relationship between in situ NH₄⁺ concentrations and chl a was assessed, both on a seasonal basis and across the entire dataset, to determine whether chl a could serve as a reliable proxy for NH₄⁺ estimation. Given the weak correlation observed between NH₄⁺ and chl a, the feasibility of direct NH₄⁺ retrieval was further evaluated. A band-selection procedure was applied to retain Sentinel-2 MSI spectral bands with the highest correlation with NH₄⁺ while minimizing multicollinearity arises from spectral “overlap” among bands. The Extreme Gradient Boosting—Particle Swarm Optimization (XGBoost-PSO), Random Forest—Bayesian Optimization with Tree-structured Parzen Estimator (RF-BO-TPE) and Support Vector Regression—Bayesian Optimization with Tree-structured Parzen Estimator (SVR-BO-TPE) models were developed for NH₄⁺ estimation. Models’ generalization ability was subsequently assessed using the difference in the coefficient of determination (ΔR²) and overfitting ratios derived from RMSE and MAE (OR_RMSE and OR_MAE).

2. Materials and Methods

2.1. Study Area

The marine cage aquaculture site is located in the southeastern part of Pagasitikos Gulf (Eastern Mediterranean, Greece), within a semi-enclosed cove (39◦07′38” N, 23◦09′ 17” E). It comprises 18 polyethylene cages, including 12 cages with dimensions of 7.5 × 15 × 8 m and 6 cages with dimensions of 15 × 15 × 10 m (length × width × depth). Sea bream (Sparus aurata) and sea bass (Dicentrarchus labrax) are intensively cultivated, using pelleted and extruded diets. The average annual standing stocks is 130 tonnes, while the final stocking density is calculated at 8 kg/m³. The feed conversion ratio (FCR), which represents the feed requirement per kilogram of body weight gain, is 1.7. The recorded circulation patterns indicate that the cove is characterized by weak currents throughout the year [19]. The prevailing circulation is directed from the northwest toward the inner part of the cove, where the fish cages are located, before exiting toward the southwest, with an average current speed of 3.5 cm s^-1 [20]. The study area does not exhibit significant anthropogenic influence, with no identifiable pollution sources detected, while the surrounding landscape is dominated by olive cultivation and shrub vegetation. Consistent with the Eastern Mediterranean pattern of primary production [21], the cove represents an oligotrophic, phosphorus (P)-limited system, where low P availability constrains phytoplankton growth and maintains chl a at low levels [19].

2.2. In Situ Measurements

The sampling was carried out over six seasons, between summer 2021 and spring 2024, using a predefined sampling grid [20]. A total of 152 surface (0m) water samples were collected by means of a 1.4 ℓ Lymnos bottle and placed in 500 ml plastic vials for NH₄⁺ estimation. Mercury chloride (2 ml of HgC_l2 0.1 M per 500 ml sample) was added to inhibit any biological activity. Then, samples were stored at -20 °C until the analysis took place. Chemical analysis of the NH₄⁺ was based on the methodology developed by [22], using a spectrophotometer (SHIMADZU UV-1800). Simultaneous surface in situ measurements of sea surface temperature (SST), sea surface salinity (SSS), pH and chl a were also obtained using a conductivity-temperature-depth (CTD) instrument (SEABIRD-19plus). All sampling procedures and measurements were carried out during the feeding period, between 09:00 and 14:00.

2.3. Statistical Analysis of In Situ Measurements

Descriptive statistics (minimum, maximum, mean ± standard deviation and median) were calculated for each WQP across seasons [23]. The Shapiro-Wilk (W) test [24] was conducted to evaluate whether the WQPs are normally distributed. The W test is widely used for assessing normality across a wide range of distributions and sample sizes [25,26] by evaluating the correlation between the ordered sample observations and the expected order statistics from a theoretical normal distribution. When this correlation is weak, the resulting p value falls below 0.05, indicating a deviation from normality [24]. The W test is given by Equation (1) [25]:

$$W={\displaystyle \frac{{\left({\sum }_{i=1}^n{a}_i{y}_i\right)}^2}{{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}}$$

(1)

where $y_i$ is the $i^{th}$ order statistic, $\overline{y}$ is the sample mean, $a_i=\left(a_1,\dots ,a_n\right)={\displaystyle \frac{m^T{V}^{-1}}{{\left(m^T{V}^{-1}{V}^{-1}m\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}}$ and $m={(m_1,\dots ,m_n)}^T$ are the expected values of the order statistics of independent and identically distributed random variables samples from the standard normal distribution and V is the covariance matrix of those order statistics. Along with the W test, histograms [23] were also created to examine the normality assumption of WQPs. The strength of association between each pair of the five WQPs was assessed using Spearman’s rank correlation coefficient (ρ). The ρ is a standardized, non-parametric statistic that can be applied when the assumption of normality is violated. It works by fist ranking the data and then applying the Pearson’s correlation coefficient (r_s) (Equation (2)) [27],

$$r_s={\displaystyle \frac{cov(R_x,R_y)}{{\sigma }_{Rx}{\sigma }_{Ry}}}$$

(2)

where $R_x,R_y$ are the ranked versions of variables x and y, cov is the covariance between ranks and σ is the standard deviation of ranks. The Spearman’s rank correlation coefficient (ρ) ranges between -1 ≤ ρ ≤ +1 and a high absolute value indicates that there is a monotonic, but not necessarily linear correlation between two variables [28]. In the present study, the correlation between WQPs was calculated both on a seasonal basis and across all seasons. In addition, scatterplots [23] were created to illustrate the relationship between each pair of variables.

The Krustal-Wallis (H) test [29] was selected to assess whether the values of each WQP differed significantly among seasons. It is a non-parametric, rank test, applied when the assumption of normality is not met, and it is used to compare three or more independent groups in terms of their distributions. The null-hypothesis of the H test assumes that the distribution of the test variable is identical across the diverse groups. Rejecting the null hypothesis means that the distributions differ in their central tendency, spread and/or variability [23]. The H test is given by the Equation (3) [29],

$$H={\displaystyle \frac{12}{N(N+1)}}{\displaystyle \sum _{i=1}^k}{\displaystyle \frac{R_i^2}{n_i}}-3(N+1)$$

(3)

where N is the total number of observations, k is the number of groups, $n_i$ is the number of observations in group i and $R_i$ is the sum of ranks in group i.

The Dunn’s post hoc test [30] was applied to identify the specific seasons between which statistically significant differences were detected for each WQP. This non-parametric pairwise multiple comparison procedure was performed after rejecting the null hypothesis of the Krustal-Wallis. The Dunn test evaluates differences between groups based on ranked observations rather than the original data values, making appropriate for datasets that do not satisfy the assumption of normality. All observations were ranked jointly across groups and the mean rank for each group was calculated by Equation (4) [31],

$${\overline{W}}_i={\displaystyle \frac{W_i}{n_i}}$$

(4)

where ${\overline{W}}_i$ is the mean rank of the ith group, $W_i$ is the sum of ranks assigned to observations within that group and $n_i$ is the number of observations in the ith group. For each pairwise comparison between two groups (e.g., A and B), the Dunn test statistics is calculated as (Equation (5)),

$$z_{AB}={\displaystyle \frac{{\overline{W}}_A-{\overline{W}}_B}{{\sigma }_{A{\rm B}}}}$$

(5)

where $z_{AB}$ is the standardized test statistic, ${\overline{W}}_A$ and ${\overline{W}}_B$ are the mean ranks of groups A and B, respectively and ${\sigma }_{A{\rm B}}$ is the standard error of the difference between mean ranks, calculated as (Equation (6)),

$${\sigma }_{AB}=\sqrt{\left[{\displaystyle \frac{N(N+1)}{12}}-{\displaystyle \frac{\sum \left({{\tau }_S}^3-{\tau }_s\right)}{12\left(N-1\right)}}\right]\left({\displaystyle \frac{1}{n_A}}+{\displaystyle \frac{1}{n_B}}\right)}$$

(6)

where N represents the total number of observations across all groups, $n_A$and $n_B$ are the sample sizes of the two compared groups and ${\tau }_s$ represents the number of observations associated with the tth tied rank. The second term inside the brackets represents the correction for tied ranks. When no tied observations occur, this correction term becomes zero. In order to evaluate WQPs distribution across different seasons, boxplots [23]were also created.

2.4. Satellite Data

Atmospherically corrected Multispectral Sentinel-2 (MSI) Level-2A (L2A) images were selected close to the sampling dates, under free cloud conditions, and provided through the Copernicus Data Space Ecosystem (CDSE) [32]. Based on the coordinates received during sampling, pixel-to-point matching was performed, and surface reflectance values were extracted using SNAP (v.13.0). A factor of 1/10.000 was applied to L2A digital numbers (DN) to retrieve physical surface reflectance values [33]. The 60 m resolution Level-2A spectral bands were excluded due to their coarse spatial resolution, since mixed pixels containing water and non-water features (e.g., land-water interfaces and floating fish cages) may reduce the reliability of spectral information used for NH₄⁺ modelling.

2.5. Spectral Band-Selection

A band-selection procedure was applied to retain Sentinel-2 MSI spectral bands with the highest correlation with NH₄⁺ while minimizing multicollinearity arising from spectral “overlap” among bands. Multicollinearity refers to a data issue that occurs when two or more predictor variables, within a dataset, are highly correlated [23]. In statistical analysis and ML, strong relationships between input variables may reduce the reliability and precision of the estimates, leading to overfitting and poor generalization to new, unseen data [34]. The procedure combined Spearman correlation-based ranking and variance inflation factor (VIF)-based spectral band selection. VIF is the most widely used indicator of multicollinearity and is obtained through a separate regression procedure [23]. VIF values ≥ 10 and 5 < VIF ≤ 10 indicate severe and moderate multicollinearity, respectively [35].

Spectral bands were first ranked based on their Spearman correlation coefficients with NH₄⁺ because the Shapiro-Wilk test indicated that the data were not normally distributed (p < 0.05). The band-selection procedure then iteratively evaluated each band, beginning with the highest-ranked one. At each iteration, the “candidate” band was temporarily added to the subset of previously selected bands and multicollinearity was assessed using the VIF. The “candidate” band was retained only when the maximum VIF of the resulting subset remained below the predefined threshold (VIF ≤ 5), otherwise, it was discarded. This procedure was repeated until all spectral bands were evaluated. To assess the effectiveness of the proposed band-selection procedure, the complete set of Sentinel-2 MSI spectral bands, at 10 m and 20 m spatial resolutions, was also used for NH₄⁺ estimation.

2.6. Machine Learning (ML) Models

The XGBoost, RF and SVR were developed for NH₄⁺ estimation. These models were selected because they demonstrate better predictive performance compared to other regression-based ML models and neural networks (NN) architectures under limited data availability [14]. The XGBoost was implemented in Python using the XGBoost library, whereas RF and SVR were implemented using the scikit-learn library.

XGBoost is an ensemble ML model designed for improved speed and performance, which use boosting and gradient descent techniques to combine multiple decision trees. In regression tasks, each new tree is built sequentially in a way that focuses on correcting the errors of the previous trees, influenced by the residuals of the earlier ones. The final prediction is obtained by aggregating the outputs of all trees [36]. RF is an ensemble ML model that applies the bagging strategy to combine multiple decision trees. In regression tasks, individual trees are trained on randomly selected subsets of the training data. The final output is computed by averaging the output of all trees, leading to improved generalization and reduced variance. SVR is a supervised ML model designed for regression tasks. The algorithm implicitly maps the input data into a higher-dimensional feature space using kernel functions, enabling the modeling of complex nonlinear relationships between input variables and the target [36,37].

Hyperparameter optimization (HPO) of the XGBoost model was performed using Particle Swarm Optimization (PSO) implemented through PySwarms library in Python, whereas the RF and SVR models were optimized using Bayesian optimization based on the Tree-structured Parzen Estimator (BO-TPE) algorithm within the Optuna framework. The hyperparameters considered for the XGBoost model included n_estimators, max_depth, learning_rate, subsample and colsample_bytree, while those considered for the RF model included n_estimators, max_depth, criterion, min_samples_split, min_samples_leaf and max_features [36]. For the SVR model, the radial basis function (RBF) was chosen as the kernel function, following previous studies that have successfully applied this kernel for similar modelling tasks [14]. Additionally, the RBF was selected due to its ability to capture complex nonlinear relationships and its relatively simple model configuration [36,37]. Therefore, the hyperparameters considered for the SVR model included C, gamma and epsilon. In all models, the optimization process was performed by minimizing the 5-fold cross-validation mean squared error (CV-MSE) [36].

2.7. Assessment of ML Models’ Generalization Ability

The dataset was randomly split into two subsets at an 8:2 ratio, with 80% of the samples (n = 121) used for training the ML models and the remaining 20% (n = 31) used for testing. The same train-test split was used for all developed models to ensure a fair and consistent evaluation of model performance. Three commonly used performance evaluation metrics, including the coefficient of determination (R²), root mean squared error (RMSE) and mean absolute error (MAE) were calculated based on the training subset. The R² provides the proportion of variance explained by the model and is described by Equation (7):

$$R^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}{{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}}$$

(7)

where n is the number of samples, $y_i$and ${\widehat{y}}_i$ represent the actual and predicted value for sample i, respectively, while $\overline{y}$ denotes the average of actual values. The R² is typically ranges between 0 and 1, whereas higher values indicate a better fit of the model. The RMSE (Equation (8)) and MAE (Equation (9)) quantify the magnitude of the models’ predictive errors and are expressed in the same units as the target variables. Models’ performance improves as the values of the error-based metrics decrease [18,38].

$$RMSE=\sqrt{{\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}$$

(8)

$$MAE={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n\left|y_i-{\widehat{y}}_i\right|$$

(9)

The developed models were tested on the unseen dataset using the same metrics. Models’ generalization ability was quantified using the difference between the R² values obtained from the training and testing datasets (ΔR²) [39] and the overfitting ratios (OR_RMSE and OR_MAE) which are given by Equations 10 and 11:

$${OR}_{RMSE}={\displaystyle \frac{{RMSE}_{test}}{{RMSE}_{train}}}$$

(10)

$${OR}_{MAE}={\displaystyle \frac{{MAE}_{test}}{{MAE}_{train}}}$$

(11)

where RMSE_test, MAE_test, RMSE_train and MAE_train are the RMSE and MAE calculated during testing and training process of ML models, respectively. ΔR² values close to 0, as well as RMSE and MAE ratios close to 1, indicate the absence of overfitting and, consequently, good generalization ability [38,40].

3. Results

3.1. In Situ Measurements

Table 1 shows the minimum, maximum and mean ± standard deviation (SD) of in situ SST, SSS, pH, chl a and NH₄⁺ measurements across different seasons. To assess the distribution of each WQP, histograms were used alongside the Shapiro-Wilk test, both on a seasonal basis and across all seasons. The histograms initially indicated deviations from the normal distribution, which were subsequently confirmed by the Shapiro-Wilk test (p < 0.05), resulting in the rejection of the null hypothesis of normality (Figure 1a–e).

The non-parametric Spearman’s rank correlation coefficient (ρ) was used to evaluate the pairwise correlation among WQPs, both seasonally and across the entire dataset. For individual seasons, the ρ values ranged from -0.26 to 0.18, indicating a weak correlation between chl a and NH₄⁺, whereas a weak negative correlation (ρ = -0.41) was observed when considering all seasons together. Scatterplots were also used to illustrate how the two variables are related to each other. In general, the data points showed a random, non-uniform, scattered pattern, without a specific direction, indicating the absence of relationship between them (Figure 1a–e).

The non-parametric Krustal-Wallis test was used to assess whether the values of each WQP differed significantly among seasons. The null-hypothesis of Krustal-Wallis test was rejected (p < 0.05), indicating significant seasonal differences for all WQPs. The Dunn’s test was then conducted to identify the specific seasons between, which statistically significant differences were detected. The results of the pairwise multiple comparisons are given in Table 2. The boxplots (Figure 2) represent WQPs distribution across different seasons.

3.2. Spectral Band-Selection

In the present study, the band-selection procedure was initiated with 10 Sentinel-2 MSI spectral bands (B02, B03, B04, B05, B06, B07, B08, B8A, B11 and B12). Firstly, the VIF was used to assess potential multicollinearity among the spectral bands. For B07, B8A, B08, B12, B06, B11 and B05 the VIF values exceeded 10, indicating severe multicollinearity, while moderate collinearity (5 < VIF ≤ 10) was observed for B03 and B02 (Figure 3). The findings highlighted the need to address multicollinearity before incorporating spectral bands into the ML models. Then, spectral bands were ranked according to their Spearman correlation with NH₄⁺ concentrations (Figure 4) and iteratively selected under a VIF ≤ 5 to address multicollinearity. Figure 5 presents the reduced dataset of spectral bands identified through band-selection procedure (B03, B04, B05 and B08) and their relative importance in NH₄⁺ estimation.

3.3. Assessment of ML Models’ Generalization Ability

The spectral bands (B03, B04, B05 and B08) selected, through band-selection procedure, were used to develop the XGBoost (PSO), RF (BO-TPE) and SVR (BO-TPE) models. To evaluate the effect of band selection, the models were also developed using the complete set of Sentinel-2 MSI spectral bands, at 10 m and 20 m spatial resolutions. Table 3 presents the generalization ability of the ML models using the difference in the coefficient of determination (ΔR²) and overfitting ratios (OR_RMSE and OR_MAE), calculating using the R², RMSE and MAE values obtained during the training and testing processes, respectively. The XGBoost (PSO) model exhibited high overfitting when using either the complete set of spectral bands or the reduced one. The ΔR² values (0.53 and 0.45, respectively) were greater than 0, while OR_RMSE (1.80) and OR_MAE (1.77) exceeded 1, indicating lower performance and higher prediction errors in the test dataset compared with the training dataset. A similar overfitting pattern was also observed for the RF (BO-TPE) model developed using the complete set of spectral bands. Only the SVR (BO-TPE) model showed moderate overfitting when trained using the complete set of spectral bands. Both the RF (BO-TPE) and SVR (BO-TPE) models demonstrated good generalization ability when developed using only the B03, B04, B05 and B08 bands (ΔR² = 0.07, OR_RMSE = 1.11, OR_MAE = 1.02 and ΔR² = 0.07, OR_RMSE = 1.10, OR_MAE = 1.01, respectively).

4. Discussion

4.1. Seasonal Dynamics of NH₄⁺ and Chl a Concentrations in an Oligotrophic Coastal Environment Influenced by Cage Aquaculture

In the present study, the relationship between in situ NH₄⁺ and chl a measurements was evaluated, on seasonal basis and across all seasons, to determine whether chl a could serve as a proxy for NH₄⁺ estimation. A weak correlation between NH₄⁺ and chl a was observed both within individual seasons and across the entire dataset.

The distinct seasonal patterns in chl a and NH₄⁺ concentrations may have contributed to the weak correlations observed. Chl a concentrations reached maximum values during winter, reflecting enhanced vertical mixing of the water column [41]. This process facilitates the upward transport of nutrients from deeper layers into the euphotic zone, thereby supporting increased phytoplankton growth and higher chl a concentrations. In contrast, lower chl a values were observed during spring and summer due to the development of water column stratification, which limits nutrient replenishment in surface waters [19]. On the other hand, NH₄⁺ concentrations reached their maximum values during summer, coinciding with the period of enhanced nutrient release associated with aquaculture operations. Nutrient release from cage aquaculture occurs continuously throughout the year, with maximum nutrient inputs generally observed during summer [42]. As water temperature increases, fish metabolic rates rise, resulting in higher feed requirements and increased feed supply, which typically peaks towards the end of August. In contrast, during winter, when water temperature decreases, fish growth is reduced and feed supply is limited to maintenance levels required to sustain the existing biomass [3].

The NH₄⁺ inputs associated with cage aquaculture did not translate into a substantial increase in chl a concentrations. This finding is consistent with the oligotrophic characteristics of the semi-enclosed cove. The limited P availability in the study area could explain the lack of a proportional increase in phytoplankton biomass, maintaining chl a at low levels, despite NH₄⁺ inputs [43,44]. In addition to P limitation, multiple factors, such as the water-column stratification and mixing, light availability, hydrodynamic processes and water renewal [44,45], as well as the rapid uptake of NH₄⁺ [46] and its oxidation through nitrification [47,48], are often influence primary production and chl a concentrations. The combined effect of these factors may have contributed to the lack of a strong relationship between chl a and NH₄⁺, illustrating the complex NH₄⁺ dynamics governing NH₄⁺ variability in the studied water body.

4.2. Effectiveness of Spearman Correlation- and VIF-Based Band Selection on Machine Learning Model Generalization and Overfitting

The weak correlation observed between NH₄⁺ and chl a indicated that chl a could not serve as reliable proxy for NH₄⁺ estimation, leading to evaluating a direct NH₄⁺ retrieval. A major challenge associated with ML models is overfitting, which occurs when models capture noise or irrelevant patterns from the training data rather than identify the actual relationships between the input and the target variables. This could result in models that perform well on the training data but fail to accurately predict new, unseen data [18].

In general, SVR is more sensitive to feature redundancy than tree-based ensemble methods such as RF and XGBoost. In SVR, additional irrelevant or noisy features may increase the risk of overfitting, causing the model to capture random fluctuations in the training data rather than the underlying relationships between the input and target variables. Moreover, redundant and correlated features may alter the geometry of the feature space and increase the complexity of the optimization problem without providing additional predictive information. Therefore, feature selection is often recommended prior to SVR training [49,50]. RF is generally robust to redundant features in terms of predictive performance. However, correlated features may bias variable importance estimates and reduce model interpretability [51,52]. In contrast, XGBoost incorporates regularization, shrinkage and feature subsampling, making it comparatively more robust to redundant predictors [53].

A band-selection procedure was applied to retain Sentinel-2 MSI spectral bands with the highest correlation with NH₄⁺ while minimizing multicollinearity arises from the spectral “overlap” among bands. The selected bands (B03, B04, B05 and B08) were used as input features to the XGBoost (PSO), RF (BO-TPE) and SVR (BO-TPE) models. To assess the effectiveness of the proposed band-selection procedure on ML models’ generalization ability, the complete set of Sentinel-2 MSI spectral bands, at 10 m and 20 m spatial resolutions, was also used for NH₄⁺ estimation. The ΔR², OR_RMSE and OR_MAE metrics were used for model evaluation. When the complete set of Sentinel-2 MSI spectral bands was used for model training, all three models exhibited moderate to high overfitting, indicating limited ability to generalize to unseen data. In contrast, when B03, B04, B05 and B08 were used as input features, two of the three models (RF-BO-TPE and SVR-BO-TPE) exhibited low overfitting, demonstrating high generalization ability and low prediction errors on the test set relative to the training set. Consequently, applying the band-selection procedure based on Spearman correlation and VIF, improved the predictive robustness of the RF and SVR models, by reducing their tendency to overfit. These findings suggest that removing redundant information may help mitigate overfitting, particularly for models that are more sensitive to feature redundancy.

4.3. Applicability and Limitations of RF (BO-TPE) and SVR (BO-TPE) Models for NH₄⁺ Estimation and Monitoring in Complex Coastal Waters Under the Influence of Cage Aquaculture

In environmental studies, careful consideration should be given to the selection of evaluation metrics used in regression-based ML modelling to ensure meaningful interpretation of model performance. Although R² is a commonly used metric that provides information on the proportion of variance explained by the model, it may not be sufficient on its own to fully capture the complex non-linear relationships between the input and the target variables [38]. In the present study, NH₄⁺ concentration variability is influenced not only by natural nutrient cycling but also by nutrient inputs from cage aquaculture, the oligotrophic and P-limited conditions of the water body and multiple additional environmental factors, which may not be fully captured by the available Sentinel-2 MSI spectral bands. Therefore, model evaluation should incorporate complementary metrics, such as RMSE and MAE, which provide information on the magnitude of prediction errors in the original measurement units. In general, the developed RF (BO-TPE) and SVR (BO-TPE) models showed limited ability to explain NH₄⁺ variability, however, the obtained RMSE and MAE values in the testing dataset (0.29 and 0.22 μM, respectively), when interpreted relative to the observed variability of NH₄⁺ concentrations (range: 0-1.46 μM, SD: ± 0.20 μM), indicated that the models provided estimates within a reasonable deviation of the actual values. Together with the observed robustness against overfitting, these findings suggest that, despite the complexity of NH₄⁺ dynamics, the developed models could provide useful estimates for NH₄⁺ monitoring, rather than complete representations of the underlying biogeochemical processes governing NH₄⁺ variability.

5. Conclusions

The use of chl a as a proxy for the retrieval of non-OACs is a common approach in eutrophic and hypertrophic water bodies. However, this study demonstrated that chl a has limited applicability for NH₄⁺ estimation, as aquaculture-driven nutrient inputs and oligotrophic, P-limited conditions resulted in a decoupling of NH₄⁺ dynamics from chl a variability. This finding highlighted the need to evaluate the potential for direct NH₄⁺ estimation. For this purpose, ML models were employed to assess the feasibility of direct NH₄⁺ retrieval, owing to their ability to capture complex non-linear relationships between spectral information and non-OACs. A band-selection procedure was applied by retaining Sentinel-2 MSI spectral bands with the highest correlation with NH₄⁺, while minimizing multicollinearity arises from spectral “overlap” among bands. The selected bands (B03, B04, B05 and B08) enhanced the robustness of the RF (BO-TPE) and SVR (BO-TPE) models by mitigating overfitting, thereby improving their ability to generalize to unseen data. This study illustrates the potential of Sentinel-2 MSI for NH₄⁺ estimation as an alternative approach to conventional monitoring methods. However, the developed ML models may require recalibration with new data before being applied to water bodies exhibiting different trophic characteristics, in order to improve their transferability.